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This  dissertation  considers  the  problem  of  source  localization  and  parameter 
estimation  with  antenna  arrays.  The  problem  is  to  estimate  the  parameters  of  the 
incident  electromagnetic  plane  waves  with  antenna  arrays.  Our  focus  is  on  array 
geometry,  sensor  characteristics,  narrowband  and  wideband  signals,  and  estimation 
algorithms  for  the  estimation  problem. 

In  particular,  the  RELAX  algorithm,  recently  proposed  for  temporal  spectral 
analysis,  is  extended  to  solve  the  spatial  problem  of  angle  and  waveform  estimation  for 
both  narrowband  and  wideband  plane  waves  arriving  at  a  uniform  linear  array.  Unlike 
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most  existing  high  resolution  algorithms,  the  narrowband  and  wideband  RELAX  al- 
gorithms are  robust  against  the  presence  of  unknown  spatially  colored  noise.  Further, 
the  wideband  RELAX  algorithm  does  not  need  the  initial  angle  estimates  that  exist- 
ing wideband  algorithms  require  to  construct  focusing  matrices,  which  causes  biases 
in  angle  estimates.  The  wideband  RELAX  algorithm  naturally  focuses  the  narrow- 
band components  in  the  spatial  frequency  domain.  Both  numerical  and  experimental 
examples  are  used  to  demonstrate  the  performance  of  the  RELAX  algorithm  and  com- 
pare the  performance  of  RELAX  with  that  of  other  well-known  algorithms  including 
ESPRIT  with  forward/backward  spatial  smoothing,  MODE/WSF,  and  AP/ANPA 
for  narrowband  signals;  and  CSM-ESPRIT  for  wideband  signals.  We  also  show  that 
better  parameter  estimates  can  be  obtained  by  using  RELAX  as  compared  to  using 
those  other  algorithms  by  means  of  both  numerical  and  experimental  results. 

To  exploit  the  advantages  of  array  geometry  and  antenna  sensors,  a  polariza- 
tion sensitive  linear  array  that  consists  of  Co-centered  Orthogonal  Loop  and  Dipole 
(COLD)  pairs  are  proposed  for  the  estimation  of  the  parameters  of  completely  polar- 
ized narrowband  electromagnetic  plane  waves.  The  performance  of  both  angle  and 
polarization  estimation  using  the  COLD  array  are  shown  to  be  greatly  improved  as 
compared  to  using  a  crossed  dipole  array.  A  MODE  algorithm  is  presented  for  both 
angle  and  polarization  estimation  of  correlated  (including  coherent)  or  uncorrelated 
incident  signals  with  a  COLD  array.  Numerical  example  are  given  to  show  the  better 
estimation  performance  of  the  MODE  algorithm  as  compared  to  those  of  the  MUSIC 
and  NSF  algothrithms.  Finally,  we  devise  a  RELAX  approach  for  angle,  polarization 
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and  waveform  estimation  of  narrowband  signals  with  a  COLD  array.  We  also  use  nu- 
merical examples  to  demonstrate  the  superior  performance  of  RELAX  as  compared 
to  MODE  when  the  additive  noise  is  spatially  colored  and  unknown. 
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CHAPTER  1 
BACKGROUND  AND  MOTIVATION 

1.1     Background 

Many  theoretical  studies  on  parameter  estimation  witli  an  array  of  sensors  have 
been  carried  out  and  deep  insight  has  been  achieved  in  the  past  two  decades.  The  re- 
search in  sensor  array  processing  was  originally  motivated  by  its  applications  in  source 
localization  and  interference  suppression  in  radar  and  sonar.  Many  algorithms  have 
appeared  in  the  literature  for  estimating  signal  parameters  from  the  measurement 
output  of  a  sensor  array. 

The  methods  in  sensor  array  processing  can  be  classified  into  two  categories: 
nonparametric  methods  and  parametric  methods.  The  nonparametric  methods  do 
not  make  any  assumption  on  the  statistical  properties  of  the  data.  Spatial  filtering 
techniques  [1,  2]  are  the  early  approaches  to  perform  a  space-time  processing  of  data 
sampled  at  an  array  of  sensors.  Beamforming  and  Capon's  methods  [3,  4,  5,  6,  7] 
are  the  typical  nonparametric  methods.  Their  idea  is  to  form  some  spectrum-like 
function  of  the  parameter (s)  of  interest,  and  then  take  the  locations  of  the  highest 
peaks  of  the  function  as  the  estimates.  All  these  approaches  have  an  inherent  limi- 
tation of  poor  resolution.  They  are  usually  used  in  situations  where  the  information 
about  the  statistical  properties  are  not  available.  However,  spatial  filtering  methods 
with  an  increasing  number  of  novel  applications  inspired  much  of  the  subsequent 


efforts  in  statistical  signal  processing.  The  well-known  Maximum  Entropy  (ME) 
spectral  estimation  method  in  geophysics  by  Burg  [8]  and  the  Yule- Walker  autore- 
gressive  estimation  method  are  the  early  parametric  approaches.  The  introduction  of 
subspace-based  techniques  [9,  10,  11]  provided  a  new  geometric  interpretation  for  the 
sensor  array  processing  problem.  The  vector  space  formulation  of  the  sensor  array 
problem  resulted  in  a  large  number  of  algorithms  [12,  13,  14,  15].  The  Maximum 
Likelihood  (ML)  parameter  estimates  can  also  be  derived  for  the  sensor  array  prob- 
lem in  an  appropriate  statistical  framework.  ML  estimation  is  a  systematic  approach 
to  many  parameter  estimation  problems,  and  has  been  studied  by  many  researchers 
[16,  17,  18,  19,  20,  21,  22].  Unfortunately,  the  ML  method  usually  requires  a  multidi- 
mensional non-linear  optimization  search  at  a  considerable  computational  complexity. 
Reduced  computational  cost  is  generally  achieved  by  the  use  of  a  suboptimal  esti- 
mator. Much  work  on  devising  a  family  of  suboptimal  estimators  and  the  analysis 
of  the  performance  of  the  estimators  has  been  done  [23,  24,  25,  26,  27,  28].  Most 
of  these  algorithms  were  motivated  by  the  subspace  based  technique  introduced  in 
the  MUSIC  algorithm.  The  subspace-based  approach  results  in  a  resolution  that  is 
not  limited  by  the  array  aperture,  provided  that  the  number  of  data  samples  or  the 
signal-to-noise  ratio  (SNR)  is  sufficient  large.  A  review  of  the  development  of  array 
signal  processing  algorithms  will  given  in  the  next  chapter. 

1.2     Motivation  and  Contributions 

Array  signal  processing  was  centered  on  the  ability  to  fuse  data  from  data  acqui- 
sition systems,  such  as  an  array  of  sensors,  to  carry  out  a  given  estimation  task.  Many 


sophisticated  estimation  algorithms  were  the  results  of  an  attempt  of  researchers  to 
go  beyond  the  classical  Fourier-limit.  However,  many  high  resolution  array  process- 
ing algorithms  are  usually  very  sensitive  to  the  presence  of  unknown  spatially  colored 
noise.  Thus  the  performance  of  those  algorithm  is  frequently  very  poor  when  used 
in  practical  applications.  To  derive  more  robust  algorithms,  we  need  to  relax  the 
additive  white  noise  assumption  made  by  most  existing  high  resolution  algorithms. 

We  propose  a  RELAX  (RELAXation  algorithm  for  the  minimization  of  the  NLS 
criterion)  algorithm  for  angle  and  waveform  estimation  of  narrowband  plane  waves 
arriving  at  a  uniform  linear  array.  The  RELAX  algorithm  is  robust  against  the 
presence  of  unknown  spatially  colored  noise.  Since  most  of  the  algorithms  in  array 
processing  are  devised  only  for  narrowband  signals,  the  algorithm  development  for 
wideband  signals  has  only  received  some  limited  attention.  Although  some  subspace- 
based  methods  [29,  30]  were  introduced  for  the  wideband  case,  most  of  them  use 
focusing  matrices,  which  may  result  in  biased  estimates.  We  propose  a  wideband 
RELAX  algorithm  for  angle  estimation  of  wideband  plane  waves  arriving  at  a  uniform 
linear  array.  Unlike  other  wideband  algorithms,  the  wideband  RELAX  does  not 
need  initial  angle  estimates  and  naturally  focuses  the  narrowband  components  in  the 
spatial  frequency  domain. 

It  is  also  important  to  take  advantage  of  array  geometries  and  receiving  proper- 
ties of  antenna  elements.  Although  many  algorithms  have  been  developed  for  array 
signal  processing  recently,  the  characteristics  of  specific  antenna  sensors  are  only 


beginning  to  attract  more  attention.  Previous  work  on  angle  and  polarization  es- 
timation using  crossed  dipoles  [31,  32]  and  orthogonal  dipoles  and  loops  [33]  are 
examples  of  using  specific  antenna  sensors  to  estimate  the  angles  and  polarizations 
of  incident  narrowband  electromagnetic  plane  waves.  To  exploit  the  advantages  of 
arra}'  geometry  and  antenna  sensors,  a  polarization  sensitive  linear  array  that  con- 
sists of  Co-centered  Orthogonal  Loop  and  Dipole  (COLD)  pairs  are  investigated  for 
the  problem  of  parameter  estimation  of  completely  polarized  electromagnetic  (EM) 
plane  waves.  The  performance  of  both  angle  and  polarization  estimation  using  the 
COLD  array  are  shown  to  be  greatly  improved  as  compared  to  using  a  crossed  dipole 
array.  We  propose  an  efficient  MODE  (Method  Of  Direction  Estimation)  algorithm 
for  both  angle  and  polarization  estimation  of  correlated  (including  coherent)  or  un- 
correlated  completely  polarized  narrowband  signals  with  a  COLD  array.  We  also 
devise  a  RELAX  algorithm  that  can  be  used  with  a  uniform  linear  COLD  array 
for  angle,  polarization  and  waveform  estimation  of  completely  polarized  narrowband 
plane  waves. 

The  high  resolution  parameter  estimation  approaches  we  develop  in  this  disser- 
tation can  be  used  for  source  localization  in  radar  and  sonar  as  long  as  incident  waves 
can  be  modeled  as  having  discrete  angles  of  arrival.  For  example,  in  radar  applica- 
tions, the  radar  returns  of  targets  impinging  on  the  antenna  array  of  a  radar  receiver 
enable  estimation  of  the  directions  of  targets  with  these  high  resolution  parameter 
estimation  algorithms.  The  parameters  of  interest  include  angles-of-arrival  (AOAs), 
waveforms,  and  polarization  states  of  the  incident  signals. 


1.3     Outline  of  the  Dissertation 

Chapter  2  provides  the  background  material  for  the  remainder  of  the  disserta- 
tion. The  material  includes  brief  reviews  of  electromagnetic  wave  propagation  and 
polarization,  and  the  sensor  array  processing  algorithms.  This  chapter  also  presents 
the  fundamental  sensor  array  data  models  for  the  dissertation. 

Chapters  describes  how  the  RELAX  algorithm  can  be  used  for  angle  and  wave- 
form estimation  of  narrowband  signals  with  a  uniform  linear  array  for  the  case  of 
multiple  snapshots.  To  evaluate  the  performance  of  the  RELAX  algorithm,  we  ap- 
ply it  to  both  simulated  and  experimental  data,  and  compare  its  performance  with 
that  of  other  well-known  algorithms  including  ESPRIT  (Estimation  of  Signal  Parame- 
ters via  Rotational  Invariance  Techniques)  with  forward/backward  spatial  smoothing, 
MODE/WSF  (Weighted  Subspace  Fitting),  and  AP  (Alternating  Projection)/ANPA 
(Alternating  Notch-Periodogram  Algorithm)  for  narrowband  signals. 

In  Chapter  4  we  extend  the  RELAX  algorithm  to  the  case  of  wideband  sources 
and  also  present  numerical  and  experimental  examples  to  illustrate  the  performance 
of  the  wideband  RELAX  algorithm  and  compare  its  performance  with  that  of  CSM 
(Coherent  Signal-subspace  Method)  based  ESPRIT. 

In  Chapter  5  an  arbitrary  linear  array  that  consists  of  Co-centered  Orthogonal 
Loop  and  Dipole  (COLD)  pairs  is  proposed  for  the  problem  of  statistically  efficient  es- 
timation of  the  parameters  of  completely  polarized  electromagnetic  waves.  A  MODE 


algorithm  is  presented  for  both  angle  and  polarization  estimation  of  correlated  (in- 
cluding coherent)  or  uncorrelated  incident  signals  with  a  COLD  array.  The  perfor- 
mance of  both  angle  and  polarization  estimation  using  the  COLD  array  are  shown  to 
be  greatly  improved  as  compared  to  using  a  co-centered  crossed  dipole  (CCD)  array. 
Several  numerical  examples  are  presented  to  compare  the  MODE  algorithm  with  the 
MUSIC  algorithm  and  the  NSF  algorithm  for  both  angle  and  polarization  estimation. 
The  Cramer-Rao  bound  (CRB)  for  the  COLD  array  is  also  compared  with  that  for 
the  CCD  array. 

Chapter  6  presents  a  RELAX  approach  for  parameter  estimation  of  narrowband 
signals  arriving  at  a  uniform  linear  COLD  array.  The  statistical  performance  of 
this  RELAX  estimator  is  compared  with  that  of  the  MODE  estimator  described  in 
Chapter  5  via  numerical  examples. 

Finally,  Chapter  7  gives  the  conclusions  and  future  work. 

Parts  of  the  original  work  presented  in  this  dissertation  have  already  been  doc- 
umented in  the  publications  [34,  35,  36,  37]. 


CHAPTER  2 
WAVE  PROPERTIES  AND  SENSOR  ARRAY  PROCESSING 

2.1     Introduction 

In  this  chapter  we  provide  the  groundwork  necessary  for  the  material  to  be  de- 
veloped in  the  subsequent  chapters.  After  briefly  reviewing  the  important  properties 
of  electromagnetic  (EM)  waves  including  wave  propagation  and  wave  polarization,  we 
introduce  the  central  problem  of  the  dissertation:  the  problem  of  locating  radiating 
sources  by  using  an  array  of  passive  sensors.  Once  the  data  model  for  the  output 
signal  of  the  receiving  sensor  array  is  formed,  the  source  location  problem  is  turned 
into  a  parameter  estimation  problem.  The  estimation  problem  has  been  investigated 
under  both  nonparametric  and  parametric  approaches. 

2.2     Propagation  of  EM  Waves 

Electromagnetic  wave  propagation  is  referred  to  the  phenomena  that  a  time- 
changing  electric  field  produces  a  time- varying  magnetic  field,  which  in  turn  generates 
an  electric  field,  and  so  on  with  a  resulting  propagation  of  energy.  The  direction  of 
the  electric  field  E  and  magnetic  field  H  are  everywhere  perpendicular.  The  most 
important  and  most  fundamental  electromagnetic  waves  are  the  transverse  plane 
waves.  In  a  plane  wave  E  and  H  lie  in  a  plane.  A  wave  of  this  type  with  both  E  and 
H  transverse  to  the  direction  of  propagation  is  called  Transverse  ElectroMagnetic 
(TEM)  wave.  In  a  medium  with  spatially  constant  permeability  fi  and  permittivity 


e  and  with  no  free  charges  and  currents,  the  Maxwell's  equations  are 


VB    =    0,  (2.1) 


VE    =    0,  (2.2) 

VxE    =    -^i—,  (2.3) 

BE 

VxH    =    t^.  (2.4) 


By  combining  the  two  curl  equations  and  making  use  of  vanishing  divergences,  we 
can  find  easily  that  each  cartesian  component  of  E  and  H  satisfies  the  wave  equation: 


V^u  -  f^e^  =  0.  (2.5) 


The  wave  equation  has  the  well-known  plane  wave  solutions 

u{t,x)  =  e^^^'^''-^'\  (2.6) 

where  the  temporal  frequency  u  and  the  magnitude  of  the  wave  vector  k  are  related 
by 

A-  =  ujy/Jie,  (2.7) 

and  X  is  the  position  vector. 

Clearly,  the  phase  variation  of  the  wave  signal  u{t,x)  includes  both  temporal  and 
spatial  variations.  Spatial  far-field  receiving  conditions,  which  implies  plane  waves, 
are  assumed  throughout  the  dissertation.  The  superposition  principle  is  valid  if  more 


than  one  wave  travel  through  a  linear  medium.  Usually  these  propagating  waves 
carry  information  from  their  sources.  The  information  may  include  source  related 
signal  parameters,  such  as  angle-of- arrival  (AOA),  signal  waveform,  signal  polariza- 
tion state,  propagation  delay,  etc.  These  signal  parameters  are  very  important  in 
many  applications.  For  source  localization,  AOA,  signal  polarization  state  and  wave- 
form are  of  special  interest.  The  efficient  estimation  of  these  parameters,  which  is 
the  essence  of  sensor  array  signal  processing,  will  be  the  topic  of  our  main  interest  in 
this  dissertation. 

2.3     Polarization  of  EM  Waves 

Polarization  describes  the  orientation  of  the  electric  field  of  a  wave.  It  is  advan- 
tageous to  employ  an  array  of  diversely  polarized  antennas  since  multiple  signals  can 
be  resolved  on  the  basis  of  polarization  as  well  as  AOA.  To  specify  the  polarization 
of  the  waves,  we  consider  sinusoidal  waves  of  the  same  frequency.  In  a  plane  wave 
traveling  along  the  positive  z  direction,  the  electric  field  generally  has  both  x  and  y 
components.  The  general  expression  for  the  electric  field  of  such  a  wave  is  then  given 

by 

E  =  (EiX  +  ^2e^"y)e-^''%  (2.8) 

where  E\  and  £'2  are  real  and  77  (  —  180°  <  //  <  -|-180°)  is  the  phase  angle  between  z 
and  y  components.  The  corresponding  magnetic  field  is 


H  =  J-i-E^e^'^x  +  Ery)e~'''\  (2.9) 
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The  phase  and  relative  ampHtudes  Ei  and  E2  determine  the  state  of  polarization.  In 
the  most  general  case  of  elliptical  polarization  the  polarization  ellipse  described  by  Ei 
and  E2e^^,  as  time  progresses,  may  have  any  orientation  as  shown  in  Fig.  2.1  .  The 
line  segment  OA  is  the  semimajor  axis,  and  the  line  segment  OB  is  the  semiminor 
axis.  The  axial  ratio  is       .  ;       ^    t 


AR 


OA 
OB 


,     {I  <AR<  00). 


(2.10) 


Figure  2.1.  Polarization  ellipse. 

We  define  /?  to  be  the  tilt  angle  of  the  ellipse  (0  <  /3  <  180°),  and  a  to  be  the 
ellipticity  angle,  which  is  given  by 


a  =  tan~\TAR),      (-45°  <  a  < +45°). 


(2.11) 
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The  ellipticity  angle  a  is  negative  for  right-handed  and  positive  for  left-handed  po- 
larization.  For  the  case  shown  in  Fig.   2.1a  is  positive.  The  parameter  7  is  defined 

as 

7  =  tan-^(|^),      (0°<7<90°).  (2.12) 

-C/l 

The  geometric  relation  of  a,  /3  and  7  to  the  polarization  ellipse  is  illustrated  in  Fig. 
2.1.  The  trigonometric  interrelations  of  a,P,'y  and  r/  are  given  by  [38,  39] 

cos  27    =    cos  2a  cos  2/9,  (2.13) 

tan  2a  ,^  -,  ,\ 

tan  7/    =     ^-^,  (2.14) 

smzp 

tan 2/?    =     tan 27  COST/,  (2.15) 

sin  2a    —    sin  27  sin?/.  (2.16) 

The  Poincare  sphere  representation  of  wave  polarization  [40]  in  Fig.  2.2  clearly  shows 
the  relationship  among  the  four  angular  variables  a,  ^3, 7  and  7/.  The  polarization  state 
is  described  by  a  point  on  a  sphere  where  the  longitude  and  latitude  of  the  point  are 
related  to  parameters  of  the  polarization  ellipse  as: 

Longitude    =    2/3,  (2.17) 

Latitude    =    2a.  (2.18) 

The  polarization  state  described  by  a  point  on  a  sphere  can  also  be  expressed 
in  terms  of  the  angle  subtended  by  the  great  circle  drawn  from  a  reference  point  on 
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Polarization  state 
M  (  a,  P )  or  (y,  t|) 


Figure  2.2.  Poincare  sphere, 
the  equator  and  angle  between  the  great  circle  and  the  equator  as: 


Great-circle  angle    =    27, 


Equator-to-great-circle  angle    =    t}. 


(2.19) 
(2.20) 


Thus,  it  is  convenient  to  describe  the  polarization  state  by  either  of  the  two  sets 
of  angles  (a,  ,5)  or  (7,?;).  The  case  when  a  =  0  corresponds  to  linear  polarization. 
The  case  when  a  =  ±45°  corresponds  to  circular  polarization,  with  left  circular 
polarization  (a  =  4-45°)  at  the  upper  pole.  One  octant  of  the  Poincare  sphere  and 
polarization  states  at  specific  points  are  shown  in  Fig.  2.3.  In  the  general  case  any 
point  on  the  upper  hemisphere  describes  a  left  elliptically  polarized  wave  ranging 
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from  pure  left  circular  at  the  pole  to  linear  at  the  equator.  Likewise,  any  point  on 
the  lower  hemisphere  describes  a  right  elliptically  polarized  wave  ranging  from  pure 
right  circular  at  the  pole  to  linear  at  the  equator. 


Left  circular  polarization 

a  =45" 


Left  elliptical  polarization 

a  =  22.5",  p  =45" 


Linear  polarization 
a  =  0"  ,    (3  =45" 


Linear  polarization 

o  o 

a=0  ,    p  =  0 


Linear  polarization 
a=0",      p  =  22.5" 


Figure  2.3.  One  octant  of  Poincare  sphere  with  polarization  states. 


EM  waves  may  be  either  completely  polarized  or  partially  polarized,  with  the 
former  being  the  most  common.  A  completely  polarized  EM  wave  is  a  special  case  of 
a  more  general  type  of  EM  wave,  i.e.,  a  partially  polarized  EM  wave.  In  other  words, 
the  polarization  state  of  a  partially  polarized  EM  wave  is  a  function  of  time  while  a 
completely  polarized  wave  has  a  fixed  state  of  polarization.  In  practical  applications 
such  as  radar  and  ionospheric  radio  [41,  42],  the  state  of  polarization  of  a  returning 
wave  received  by  a  radar  with  polarization  diversity  can  vary  even  though  the  original 
transmitted  wave  is  completely  polarized. 
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2.4     Narrowband  Signals 

2.4.1     Array  Data  Model  and  Problem  Formulation 

A  general  two-dimensional  isotropic  sensor  array  system  is  shown  in  Figure  2.4. 
The  wave  field  of  the  sources  travels  through  space  and  is  sampled,  in  both  space 
and  time,  by  the  sensor  array.  Assume  that  the  array  is  planar,  each  of  which  has 
coordinate  r;  =  (x/,?/;)  and  an  impulse  response 


hi{t,vi)  =  ai{e)8it)6iri),      /=1,---,L. 


(2.21) 


Y! 


Source  1 


Sensor  Array 


Source  2 


ff 


X 


Figure  2.4.  An  arbitrary  two-dimensional  array. 


For  K  signals  impinging  on  an  array  of  L  sensors,  we  can  define  an  L  x  K 
impulse  matrix  H(^,^)  from  the  impinging  emitter  signals  with  parameter  0  = 
[  9i     $2     •••   Gk  1'  ^°  ^^^  sensor  outputs.    The  Ikth  element  of  H(^,0)  takes  the 
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form 

H/,(^, 6)  =  ai{ek)6it)8{ri),     I  =  1,- ■■  ,L,k  =  1,2,- ■■  J<.  (2.22) 

From  a  convolution  operation,  the  sensor  outputs  can  be  written  as 


K 

I 

fc=i 


y/(0  =  E«K^fc)e^'^"*-''^"^   *    .  (2-23) 


for  purely  exponential  signals  and 


K 

I 


yi{t)  =  Y:  ai{e,)e^^'^'-''^''hk{t)  (2.24) 


for  narrowband  signals,  where  u  is  the  common  center  frequency  of  the  signals  and 
Sk{t)  represents  the  complex  amplitude  (waveform)  of  the  ^th  signal.  The  carrier  e^'^'^ 
is  usually  removed  from  the  sensor  output  before  sampling. 

It  is  clear  that  the  geometry  of  a  given  array  determines  the  relative  delays 
of  the  various  angles-of-arrival  (AOA).  The  formulation  can  be  straightforwardly 
extended  to  arrays  where  additional  dimensions  provide  the  flexibility  for  more  signal 
parameters  per  source,  such  as  a  polarization  sensitive  sensor  array  for  both  AOA 
and  polarization  estimation.  We  will  introduce  a  polarization  sensitive  sensor  array 
later  in  this  dissertation. 

The  received  signal  plus  noise  gives  the  outputs  of  the  sensor  array  in  the  form 
of 


yW  = 


a(^i),    a(^2),    •••,    ^iOx) 


s{t)  +  n{t)  =  As{t)  +  n{t),  (2.25) 


where 
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a(^fc)  = 


aa(^fc)e^k^S    G2(^fc)e^■k^^     .••,    a^l^^e^'^^"'^ 


(2.26) 


<t) 


Si{t),      S2{t),---,      SK{t) 


(2.27) 


and  n(^)  is  the  noise  vector.  The  vector  a{Ok)  is  I'eferred  to  as  the  ^th  array  propa- 
gation steering  vector.  It  describes  via  k^r;  how  a  plane  wave  impinges  at  the  array 
(i.e.  AOA)  and  via  ai{9k)  how  the  sensors  affect  the  signal  ampUtude  and  phase. 
For  polarization  sensitive  sensors,  the  ai[9k)  is  related  to  the  polarization  state  of 
the  incident  electromagnetic  plane  wave.  For  omnidirectional  isotropic  sensors,  the 
ai(9k)  is  a  constant.  We  consider  both  isotropic  sensor  arrays  (Chapters  3  and  4) 
and  polarization  sensitive  sensor  arrays  (Chapters  5  and  6)  in  this  dissertation.  The 
reviews  of  previous  work  presented  below  are  for  an  isotropic  sensor  array  unless  we 
point  out  a  polarization  sensitive  array. 

For  a  uniform  linear  array  (ULA)  shown  in  Figure  2.5,  the  steering  vector  a(^;i:) 
has  the  form 


eiiOk)  = 


-\T 


1        pj^k 


j(L-i)nfc 


(2.28) 


where 


u> 


Ofc  =  —6sm(6k) 
c 


(2.29) 


is  called  the  spatial  frequency. 
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Source 


Ssine 


Sensor 


X 


Figure  2.5.  A  uniform  linear  array. 

We  can  see  from  (2.29)  that  the  vector  a{0k)  is  uniquely  defined  if  and  only  if 
Clk  is  constrained  as  \^k\  <  i".  The  condition  is  satisfied  if 


S< 


A 


(2.30) 


The  collection  of  these  steering  vector  over  the  parameter  space  of  interest 


A^{aiek)\ekee} 


(2.31) 


is  often  called  the  array  manifold.  The  parameterization  of  A  is  assumed  known.  Let 
us  define  the  set  A'^  as  the  collection  of  all  distinct  array  manifold  vectors 


^«  =  {A|A  =  [  a(^i)    ...    a(^K-)],^i<^2<---<^K}.  (2.32) 
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Hence,  A^""  is  parameterized  by  the  parameter  vector  0  =  [  0^^  •••,  0j^-  ]'^ ■  The 
array  is  assumed  to  be  unambiguous.  In  other  words,  any  A  e  A^  has  full  rank. 

The  sensor  outputs  are  appropriately  sampled  at  ^  =  1,2,  •  •  • ,  A'^  time  instances 
and  these  snapshots  y(l),  y(2),  •  •  •,  y(A^)  can  be  viewed  as  a  multichannel  ran- 
dom process,  which  is  assumed  Gaussian  in  this  dissertation.  The  characteristics  of 
Gaussian  processes  can  be  well  understood  from  its  first  and  second  order  statistics 
determined  by  the  underlying  signals  as  well  as  noise.  The  problem  of  central  interest 
for  source  localization  is  to  estimate  the  AOAs,  waveforms  (and  polarization  states 
if  a  polarization  sensitive  array  is  used)  of  emitter  signals  impinging  on  a  receiving 
array  when  a  set  of  sample  data  {y(l),y(2),  •  •  •  ,y(A'^)}  is  given. 

We  first  make  some  assumptions  for  the  additive  noise  n(^)  and  signal  waveforms 
s(^).  The  noise  vector  n(^)  is  assumed  to  be  a  stationary,  temporally  white,  zero-mean 
and  circularly  symmetric  with  unknown  covariance  matrix  Q: 

E{n{h)n"{t2)}^QSt,,t„  (2.33) 

E{n{h)n^{t2)}  =  0,  (2.34) 

where  St^^t2  is  the  Kronecker  delta,  (•)^  represents  the  complex  conjugate  transpose, 
and  {•)^  denotes  the  transpose.  Note  that  the  problem  of  angle  estimation  is  ill- 
defined  for  an  arbitrary  noise  field  without  knowledge  of  the  signal  waveform.  The 
eigenstructure-based  estimation  methods  assume  the  case  of  spatially  white  noise, 
i.e.,  Q  =  aH. 
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The  signal  waveform  s(^)  is  assumed  to  be  deterministic  unknown  or  random. 
For  the  latter  case,  the  source  covariance  matrix  is  defined  as 

R,  =  E{s(Os^(i)}.  (2.35) 


For  the  former  case, 

Rs^^Es{t)s"{t).  (2.36) 

-'^  t=i 

For  the  case  of  spatially  white  noise,  the  array  covariance  matrix  has  the  form 


R  =  E{y{t)y"it)}  =  AR,A"  +  a^I.  (2.37) 

The  eigendecomposition  of  R  results  in  the  representation 

K=J2  ^'-e.-ef  =  E,A,Ef  +  E„A„E^,  (2.38) 

where  Ai  >  •  •  •  >  A^  >  Aft-+i  =  •  •  •  =  A^,  =  cr^.  The  matrix  E^  =  [  e^,  •  •  • ,  e^-  ] 
contains  the  K  eigenvectors  corresponding  to  the  largest  eigenvalues.  The  range  space 
of  Es  is  called  the  signal  subspace.  Its  orthogonal  complement  is  the  noise  subspace 
and  is  spanned  by  the  columns  of  E„  =  [  e^'+i,    •  •  • ,    e£,  ]• 

The  eigendecomposion  of  the  sample  covariance  matrix  R  is  given  by 


R  =  T7  E  ylOy^'C^)  =  E.A.Ef  +  E„A„E^.  (2.39) 

-''    4=1 
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The  number  of  signals  K  is  assumed  known  throughout  the  dissertation.     Signal 
enumeration  methods  can  be  found  in,  e.g.,  [43,  44,  45,  46,  47,  48]. 

In  the  next  two  subsections,  we  introduce  the  most  well  known  estimation  tech- 
niques, classified  as  nonparametric  and  parametric  methods. 

2.4.2     Nonparametric  Methods 

The  nonparametric  methods  do  not  make  any  assumption  on  the  statistical 
properties  of  the  data.  The  fundamental  idea  of  the  nonparametric  methods  is  to 
form  some  spectrum-like  function  of  the  parameter(s)  of  interest,  and  then  take  the 
location  of  the  highest  peaks  of  the  function  as  the  estimates. 

The  beamforming  is  to  "steer"  the  array  in  one  direction  at  a  time  and  measure 
the  array  output  power.  The  locations  of  maximum  power  yield  the  AOA  estimates. 
The  array  is  steered  by  forming  a  linear  combination  of  the  sensor  outputs 

yF{t)  =  Y.^o*y.(i)  =  ^"yi^)-  (2-40) 

!  =  1 

Given  {yF(0}ili'  ^^^  output  power  is  measured  by 

^  E  bHO I'  =  ^  E  w^y(Oy''(^)w  =  w^Rw.  (2.41) 

Different  choices  of  weighting  vector  w  leads  to  different  beamforming  approaches. 
There  are  two  types  of  beamformers. 
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Conventional  Beamformer 

The  weighting  vector  for  conventional  beamformer  is  given  by 

WBF  =  a(^),  (2.42) 

which  can  be  interpreted  as  a  spatial  filter  matched  to  the  impinging  signal. 

Substituting  the  weighting  vector  (2.42)  into  (2.41),  the  beamforming  AOA 
estimates  are  given  by  the  locations  of  the  largest  peaks  of  the  spatial  spectrum 

Pbf{0  =  a^(6')Ra(^).  (2.43) 

Clearly,  for  a  ULA,  the  beamforming  is  a  spatial  extension  of  the  classical  peri- 
odogram  in  temporal  time  series  analysis  [7].  Thus,  the  spatial  spectrum  suffers 
from  the  same  resolution  problem  as  does  the  periodogram.  This  sets  a  limit  on  the 
resolution  achievable  by  beamforming,  which  is  27r//y  for  a  ULA. 

Capon's  Beamformer 

A  well-known  method  was  proposed  by  Capon  [3],  whose  beamformer  attempts 
to  minimize  the  power  received  from  noise  and  any  signals  coming  from  other  direc- 
tions than  6,  while  maintaining  a  fixed  gain  in  the  direction  6.  The  Capon  spatial 
filter  design  problem  is  posed  as 

minw^Rw     subject  to     w"a{9)  =  1.  (2.44) 
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The  optimal  w  is  given  by 


wcp  =  ^ ,  (■'•'±5) 


which  leads  to  the  following  spatial  spectrum  upon  insertion  into  (2.41) 


CP 


a^(^)R-ia(^) 


(2.46) 


Thus,  the  Capon  AOA  estimates  are  obtained  by  the  locations  of  the  largest  peaks 
of  the  spatial  spectrum  given  by  (2.46). 

Capon's  beamformer  outperforms  the  conventional  one  because  the  former  uses 
every  single  degree  of  freedom  to  concentrate  the  received  energy  along  the  direc- 
tion of  interest.  However,  the  resolution  capability  of  the  capon  beamformer  is  still 
dependent  on  the  array  aperture.  Since  nonparametric  methods  do  not  assume  any- 
thing about  the  statistical  properties  of  the  data,  they  can  be  used  in  situations 
where  we  lack  information  about  these  properties.  Several  alternative  methods  for 
beamforming  have  been  proposed  for  addressing  various  issues  [4,  5,  6]. 

2.4.3     Parametric  Methods 

When  some  of  the  statistical  information  of  the  data  is  available,  the  use  of 
a  nonparametric  approach  is  often  associated  with  a  degradation  in  performance 
as  compared  to  the  use  of  model-based  parametric  approach.  The  most  important 
model-based  approach  is  the  Maximum  Likelihood  (ML)  method.  The  ML  estimates 
are  those  values  of  the  unknown  parameters  that  maximize  the  likelihood  function. 
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There  are  two  different  ML  approaches  in  the  sensor  array  problem,  depending  on 
the  model  assumption  about  the  signal  waveform.  The  unknown,  deterministic  signal 
model  leads  to  the  deterministic  ML  method  [18,  17,  49,  50,  51],  while  the  Gaus- 
sian random  signal  model  results  in  the  stochastic  ML  method  [19,  52,  53].  Al- 
though the  ML  approaches  are  optimal  to  the  sensor  array  problem,  these  techniques 
are  often  deemed  exceedingly  complex.  For  the  case  of  spatially  white  noise,  the 
subspace-based  techniques  provide  many  high  resolution  and  computationally  effi- 
cient algorithms.  These  algorithms  include  MUSIC  [10,  11]  and  its  modified  versions 
[54,  55,  56,  57],  and  ESPRIT  [58,  59]  as  well  as  root-MUSIC  [60,  61]  for  a  uniform 
linear  array.  The  subspace-based  techniques  also  lead  to  the  subspace-based  ML  ap- 
proximation approaches,  which  include  Signal  Subspace  Fitting  (SSF)  [23,  62,  63] 
and  related,  yet  more  efficient  Method  Of  Direction  Estimation  (MODE)  [25,  26],  as 
well  as  Noise  Subspace  Fitting  (NSF)  [25,  64,  65].  Compared  to  the  exact  ML  ap- 
proach, the  subspace-based  ML  approximation  approaches  are  computationally  more 
attractive.  However,  the  subspace-based  ML  approximation  approaches  are  usually 
sensitive  to  the  violation  of  white  Gaussian  noise  assumption,  and  the  number  of 
incident  signals.  In  the  next  three  subsections  we  will  briefly  review  these  important 
parametric  methods. 

ML  Methods 

The  ML  methods  include  the  deterministic  ML  and  stochastic  ML. 
•  Deterministic  ML 
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The  deterministic  ML  method  maximizes  the  conditional  likelihood  function 
(given  the  signal  waveforms  and  the  incident  angles)  of  the  data.  It  is  also  called 
the  nonlinear  least-squares  (NLS)  method.  This  leads  to  the  minimization  of  the 
following  function  [18,  17]: 


NLF  =  E[y(n)  -  A{0)s{n)f[y{n)  -  A{0)s{n)].  (2.47) 

n=l 


Minimizing  the  cost  function  with  respect  to  s(l),  s(2),  •  •  •,  s(N)  gives 


s(0  =  (A^A)-iA^y(i),     t  =  l,2,---,N.  (2.48) 


By  substituting  (2.48)  in  (2.47),  we  get 


^    =    ^j:\\{I-MA"A)-'A"}y{t)r  (2.49) 

=    argmaxTr{A(A^A)-iA^R},  (2.50) 

0 


where  R  is  the  sample  covariance  matrix. 

Note  that  a  nonhnear  /i^-dimensional  optimization  problem  must  be  solved  for 
the  DML  estimator.  A  good  initial  estimate  is  important  to  guarantee  the  desired 
global  minimum  since  the  criterion  function  often  possesses  a  large  number  of  local 
minima.  A  spectral-based  method  is  a  natural  choice  for  an  initial  estimator  provided 
all  sources  can  be  resolved.  The  alternating  projection  (AP)  [66]  and  alternating 
notch-periodogram  algorithm  (ANPA)  [67]  are  two  similar  efficient  approaches  to 
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minimize  the  nonlinear  least-squares  criterion.    The  IQML  algorithm  [49,  50,  51] 
is  an  iterative  procedure  to  minimize  the  deterministic  ML  criterion  in  (2.50)  for 
a  uniform  linear  array.    The  RELAX  algorithm  presented  in  the  next  chapter  is  a 
FFT-based  method  to  minimize  the  NLS  criterion  in  (2.49). 
•  Stochastic  ML 

The  Gaussian  signal  model  assumption  not  only  is  a  way  to  obtain  a  tractable 
ML  method,  but  also  is  often  motivated  by  the  Central  limit  Theorem.  Assume  the 
signal  waveforms  are  zero-mean  with  second-order  properties 

E{s(ii)s^(i2)}  =  RA,t2,      E{s(ii)s^(i2)}  =  0.  (2.51) 

Under  the  assumption,  the  observation  process  y{t)  is  a  white,  stationary,  zero-mean 
and  circularly  symmetric  Gaussian  random  vector  with  covariance  matrix 

R  =  A{e)RsA"{B)  +  a^I.  (2.52) 

The  negative  log-likelihood  function  (ignoring  constant  terms)  of  the  complete 
data  set  y(l),  •  •  • ,  y(A'^)  is  given  by 

l{e,  R„  a^)  =  log|R|  +  tr{R-iR}.  (2.53) 

Minimization  with  respect  to  cr^  and  R^  leads  to  the  following  expressions: 

RsiB)  =  At(R  -  a-^(^)I)At^,  :      (2.54) 
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a\e)  =  -^tr{PiR}.  (2.55) 


where 

At  =  (A^A)~'a^,  (2.56) 

Pi  =  I  -  A  (a^ a)  ''  A".  (2.57) 

Substituting  (2.54)  and  (2.55)  into  (2.53)  results  in  the  following  compact  form 

of  the  criterion: 

0  =  arg  min  {log  | AR,(6>)A^  +  a^O)!]]  .  (2.58) 

0 

The  dimension  of  the  parameter  space  is  reduced  substantially  for  the  SML 
estimator.  However,  the  criterion  function  in  (2.58)  is  a  highly  non-linear  function  of 
0  ,  and  the  minimization  0  cannot  be  found  analytically  in  general.  A  Newton-type 
implementation  of  numerical  search  may  result  in  a  very  good  accuracy  when  the 
global  minimum  is  attained.  In  fact  the  stochastic  ML  estimator  has  been  shown 
to  have  a  better  large  sample  accuracy  than  the  corresponding  deterministic  ML 
estimator  [20,  62]  regardless  of  the  actual  distribution  of  the  signal  waveforms. 

Subspace-Based  Methods 

The  tremendous  interest  in  the  subspace  approach  is  mainly  the  result  of  the  in- 
troduction of  the  Multiple  Signal  Classification  (MUSIC)  algorithm  [68,  60],  although 
Pisarenkos  work  (a  special  case  of  MUSIC)  [69]  in  time  series  analysis  was  pubHshed 
in  early  70's.  These  methods  assume  that  the  additive  noise  is  either  spatially  white 
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or  its  covariance  matrix  is  known.   Our  discussions  on  the  subspace-based  methods 
below  assume  that  the  noise  is  spatially  white. 
•  MUSIC 

The  MUSIC  algorithm  is  based  on  the  fact  that  the  noise  eigenvectors  in  E„  are 
orthogonal  to  A,  i.e., 

E^a(^)  =  o,    ee{e„---,eK}.  (2.59) 

For  unique  AOA  estimates,  the  array  is  usually  assumed  to  be  unambiguous.  In 
particular,  a  ULA  is  unambiguous  if  —90°  <  ^  <  90°  and  d  <  A/2.  In  practice,  the 
MUSIC  "spatial  spectrum"  is  defined  as 


P„(,)  =       -^"W^W      .  (2.60) 


The  MUSIC  AOA  estimates  are  given  by  the  locations  of  the  largest  peaks  of 
the  "spatial  spectrum"  in  (2.60)  provided  E„E^  is  close  to  E„E^. 

In  contrast  to  the  beamforming  techniques,  the  MUSIC  algorithm  provides  es- 
timates of  an  arbitrary  accuracy  if  the  data  samples  are  sufficiently  large  or  the  SNR 
(signal-to-noise  ratio)  is  adequately  high.  A  number  of  modifications  of  MUSIC  have 
been  proposed  to  improve  or  overcome  some  of  its  shortcoming  in  various  specific 
scenarios.  The  Min-Norm  algorithm  [54,  55,  56,  57],  a  weighted  MUSIC,  exhibits  a 
better  resolution  than  the  original  MUSIC  algorithm  [70]. 
•  ESPRIT 
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A  uniform  linear  array  steering  matrix  has  the  structure  of  a  Vandermonde 
matrix  as  follow 


A  = 


J^i 


J^K 


J(i-l)fil 


J{L-1)Qk 


(2.61) 


Making  use  of  this  special  structure  results  in  several  computationally  and  statis- 
tically efficient  subspace-based  algorithms.  This  type  of  algorithm  may  include 
ESPRIT  (Estimation  of  Signal  Parameters  via  Rotational  Invariance  Techniques) 
[58,  59],  and  Root-MUSIC  [60]  that  will  be  introduced  later. 

The  ESPRIT  algorithm  exploits  a  so-called  shift  structure  of  the  Vandermonde 
matrix  A.  let  Ai  and  A2  be  the  sub-matrices  by  deleting  the  first  and  last  rows  from 
A  respectively.  Then,  Ai  and  A2  are  related  by 


Ai*, 


(2.62) 


where 


$  =  diag{e-'"\eJ"%---,e^"^-} 


J^K  ■ 


(2.63) 


Consider  the  structure  of  the  eigendecomposition  of  the  array  covariance  matrix  R 
given  in  (2.38).  If  the  signal  covariance  matrix  R^  has  rank  K  ,  then  the  matrix  E., 
will  span  a  K  -dimensional  subspace  of  A.  This  observation  implies  that  there  exists 
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a  full-rank  K  x  K  matrix  such  that 


E,  =  AT.  (2.64) 


Deleting  the  first  and  last  rows  of  (2.64),  respectively,  gives 


E,i  =  AiT,     E,2  =  A2T.  (2.65) 

Combining  (2.62)  and  (2.65)  leads  to 

E,2  =  E,i*,  (2.66) 

where  ^  =  T~^$T.  Clearly,  '9  and  $  have  the  same  eigenvalues.  The  eigenvalues 
are  given  by  e-'^* ,  A;  =  1, 2,  •  •  • ,  A' ,  which  are  related  to  AOAs.  There  are  two  different 
ESPRIT  algorithms,  depending  on  how  to  approximate  the  relation 

E«2  =  E,i*.  (2.67) 

Solving  the  approximation  relation  (2.67)  in  a  least-squares  sense  results  in  LS- 
ESPRIT,  while  in  a  Total-Least-Squares  sense  leads  to  TLS-ESPRIT. 
•  Root-MUSIC 
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The  idea  of  the  Root-MUSIC  dates  back  to  Pisarenko's  method  [69].    If  the 
signals  covariance  matrix  has  full  rank,  then  the  polynomial 

p,(z)  =  ef^a(z),      l  =  K  +  h---,L,  (2.68) 

where  a{z)  =  [1,  5,  •  •  •  ,2^"^]^,  has  K  zero-roots  at  e-'"*,  k  =  1,2,  •••,/<'.  Consider 
the  zeros  of  the  MUSIC  function 

||E^a(2)||^==a^(.)E„E^a(2).  (2.69) 

The  search  for  zeros  is  complicated  since  the  function  is  not  a  polynomial  in  2.  In 
fact,  we  are  only  interested  in  values  of  z  on  the  unit  circle.  This  suggests  that  we 
can  use  a'^{z~^)  to  substitute  for  a^(z),  which  gives  the  Root-MUSIC  polynomial 

p{z)  =  z^a'{z-')Er.E^R{z).  (2.70) 

The  polynomial  p{z)  has  2/i  roots  with  A'  roots  inside  the  unit  circle  and  K  mirrored 
images  outside  the  unit  circle.  Among  those  inside,  the  phases  of  the  K  closest  to 
the  unit  circle  gives  the  AOA  estimates. 

The  Root-MUSIC  has  empirically  been  found  to  have  a  significant  better  per- 
formance than  MUSIC  in  the  small  sample  case  according  to  [28,  71]. 
•  Coherent  Signals 


31 

For  highly  correlated  or  coherent  signals,  a  rank  deficiency  occurs  in  the  source 
covariance  matrix  Rj.  This  results  in  a  divergence  of  a  signal  eigenvector  into  the 
noise  subspace.  Thus  the  property  (2.59)  no  longer  holds  and  the  subspace  method 
fails  to  yield  consistent  estimates.  The  forward-backward  (FB)  averaging  and  spatial 
smoothing  techniques  can  be  used  with  subspace  methods  for  ULAs  in  the  limiting 
case  of  coherent  signals  [6,  72,  44,  73].  The  idea  of  the  FB  averaging  is  to  form 
a  so-called  backward  array  covariance  matrix  R^  and  then  average  the  usual  array 
covariance  R  and  R^-  The  backward  array  covariance  matrix  has  the  form 

Rb  =  JR*J  =  A^-^^-^'R.^-^-^-'^A^  +  a%  (2.71) 

where  J  is  an  L  x  1/  permutation  matrix,  whose  components  are  zero  excepts  for  ones 
on  the  anti-diagonal.  Then  the  FB  array  covariance  matrix  is  obtained  by 


Rfb  =  -(R  +  Rb)  =  AR,A^  +  a%  (2.72) 


where  R,  =  (R,  +  *-(^-^^R,$-(^-^^)/2  has  full  rank. 

The  idea  of  the  spatial  smoothing  technique  is  to  split  the  ULA  into  a  number 
of  overlapping  subarrays.  The  outputs  of  the  subarrays  can  therefore  be  averaged 
for  computing  the  covariance  matrix.  Since  the  spatial  smoothing  induces  a  random 
phase  modulation  similar  to  (2.72),  the  signals  are  decorrelated  and  results  in  a  full 
rank  covariance  matrix. 
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Subspace-Based  ML  Approximations 

Subspace-based  methods  offer  significant  performance  improvement  in  compar- 
ison to  conventional  beamforming  methods.  Although  the  MUSIC  method  yields 
estimates  with  the  same  large  sample  accuracy  as  that  of  the  deterministic  ML 
method  provided  the  signals  are  uncorrelated  [20],  it  usually  suffer  from  a  large  bias 
for  finite  samples  and  is  unable  to  cope  with  coherent  signals.  Several  parametric 
subspace-based  methods  that  practically  have  the  same  statistical  performance  as 
the  ML  method  but  less  computational  complexity  than  ML  have  been  proposed 
[25,  26,  23,  63]. 
•  Signal  Subspace  Fitting 

Recall  the  relation  in  (2.64).  Since  6  and  T  are  unknown,  it  appears  there  are 
no  values  of  0  such  that  Eg  =  AT  when  only  an  estimate  E^  of  E^  is  available.  A 
natural  solution  is  to  find  the  best  weighted  least-squares  fit  of  the  two  subspaces  by 
the  following  criterion: 


{0,T}  =  argmin||E,  -  AT||^,  (2.73) 


where  ||A||^  =  tr{AWA^}  and  W  is  a  K  x  K  positive  definite  weighting  matrix. 
This  is  a  separable  nonlinear  least-squares  problem  [74].  Substituting  the  pseudo- 
inverse  solution,  T  =  ATEs,  into  (2.73)  leads  to  the  separated  criterion  function 


^5SF  =  argminTr{PiE,WEf}.  (2.74) 

0 
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Different  choices  of  W  lead  to  different  methods.  The  choice  W  =  I  yields  Cadzow's 
method  [75].  A  statistical  analysis  suggests  other  weightings,  such  as 


W  =  (A,  -  cr'lYA,   , 


(2.75) 


leading  to  superior  performance  [25,  23,  76].  The  estimator  defined  by  (2.74)  with 
weights  given  by  (2.75)  is  called  the  MODE/WSF(Weighted  Subspace  Fitting)  method. 
For  a  ULA,  the  MODE  algorithm  [25,  26]  can  be  implemented  efficiently.  The 
idea  is  to  reparameterize  the  projection  matrix  P^  using  a  basis  for  the  nullspace 
of  A^.  The  projector  P^  above  may  be  reparameterized  in  terms  of  the  coefficients 
b  =  [  |,jj    l)^     ...    ij^.  ]-^  of  a  polynomial  defined  as 


fC  K 


K 

I 

A:=0 


A;=l 


(2.76) 


Let  B'^  be  the  following  {L  -  K)  x  L  matrix 


B^  = 


hK 


0 


61 


hK 


bo 


h 


bo 


(2.77) 


Then  P^  =  Pb,  where  Pb  denotes  the  orthogonal  projector  onto  the  range  space  of 
B.  Thus  the  reparameterized  SSF  criterion  for  MODE  is  given  by 


b  =  argminTr{B(B''B)-'B''E,WEf  }. 

b 


(2.78) 
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Clearly,  the  resulting  non-iterative  procedure  includes: 

1)  Determine  b  by  minimizing  (2.78)  with  (B^B)~^  set  to  I; 

2)  Determine  an  asymptotically  efficient  estimate  b  as 


-1, 


b  =  argminTr{B(B^B)     B^E.WEf},  (2.79) 

b 


where  (B^B)  is  formed  by  the  initial  estimate  of  b  obtained  in  Step  1; 

3)  Compute  the  roots  of  the  estimated  polynomial. 

This  MODE  algorithm  for  a  ULA  is  attractive  because  of  the  fact  that  it  is 
computationally  and  asymptotically  efficient  ,  and  it  remains  so  for  highly  correlated 
or  coherent  signals.    Two-dimensional  spatial  frequency  and  angle  estimation  using 
MODE  has  been  presented  in  [46,  77]. 
•  Noise  Subspace  Fitting 

Inspired  by  MUSIC,  if  signal  covariance  Rj  has  full  rank,  then  the  following 
relation  holds: 

Ef  A  =  0.  (2.80) 

Given  an  estimate  of  En,  a  natural  estimate  of  B  can  be  obtained  by  minimizing 
the  following  Noise  Subspace  Fitting  (NSF)  criterion: 

^Ar5F  =  argmintr{UA^E„EfA},  (2.81) 

u 

where  U  is  a  K  x  K  positive  semi-definite  weighting  matrix.    Different  choices  of 
weighting  matrix  U  result  in  signal  parameter  estimates  with  different  asymptotic 
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properties.  It  has  been  shown  in  [76]  that  the  estimates  calculated  by  (2.81)  and 
(2.74)  are  asymptotically  equivalent  if  the  weighting  matrix  is  given  by 

U  =  A^E^WEf  At^.  (2.82) 

Compared  with  SSF,  however,  NSF  has  some  inherent  disadvantages,  namely  that 
it  may  not  cope  with  coherent  signals,  and  a  two-step  procedure  has  to  be  adopted 
because  the  weighting  matrix  depends  on  6. 

The  parametric  high  resolution  algorithms  introduced  in  this  subsection  are  very 
sensitive  to  the  presence  of  unknown  spatially  colored  noise.  The  performance  of  these 
algorithms  is  usually  very  poor  when  used  in  practical  applications.  To  derive  more 
robust  algorithms,  we  may  need  to  relax  the  spatially  white  noise  assumption  made 
by  these  algorithms.  A  robust  RELAX  algorithm  will  be  proposed  in  Chapter  3  for 
angle  and  waveform  estimation  of  narrowband  plane  waves  arriving  at  a  ULA.  The 
RELAX  algorithm  does  not  assume  that  the  noise  is  spatially  white. 

2.5     Wideband  Signals 

The  methods  introduced  above  are  essentially  limited  to  narrowband  signals.  For 
wideband  signals,  the  delay  spread  is  close  to  If BW ^  where  BW  is  the  bandwidth  of 
signals.  The  wideband  signal  is  usually  decomposed  into  narrowband  components  by 
filtering  or  DFT  (Discrete  Fourier  Transform),  and  then  somehow  combine  the  differ- 
ent frequency  bins  before  using  the  narrowband  algorithms.  The  question  of  how  to 
best  combine  the  information  spread  across  different  frequencies  has  received  much 
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attention.  The  algorithms  proposed  in  [78,  79]  process  each  frequency  bin  separately 
by  using  MUSIC,  and  then  averages  the  spatial  spectrums,  whereas  [80]  proposes 
to  use  ESPRIT.  These  methods  are  incoherent  in  the  sense  that  the  directional  in- 
formation is  spread  over  a  wideband  of  frequencies,  and  each  frequency  is  processed 
separately  followed  by  some  statistical  average  of  the  results.  This  average  of  the 
results  tends  to  be  ineffective  at  low  SNR  and  thus  the  resolution  of  these  methods 
tends  to  be  limited  by  the  resolution  of  a  single  frequency  bin.  A  number  of  focusing 
algorithms  have  been  proposed  including  coherent  signal  subspace  methods  (CSM) 
[29,  30],  and  the  spatial  resampling  method  [81].  These  coherent  methods  are  an 
alternative  to  the  incoherent  methods.  These  methods  preprocess  the  data,  focusing 
the  different  frequency  components  to  a  signal  reference  frequency,  removing  the  fre- 
quency dependence  of  the  data.  Other  wideband  AOA  estimation  algorithms  have 
been  developed  in  [82,  83]. 

We  now  briefly  review  the  wideband  array  data  model  and  the  concept  of  CSM. 
Assume  that  K  incident  wideband  signals  impinge  at  a  uniform  linear  array  with  L 
sensors.  The  data  model  in  frequency  domain  has  the  form 

y(/„)  =  A(0,/„)r(/„)  +  N(/„),     n-l,2,---,J,  (2.83) 

where  y(/„)  denotes  the  L  x  1  data  vector  of  the  n-th  frequency  bin,  A(0,/„)  = 
[  af^i    f  )     •  ■  ■    a-iOw   f  )  ]  is  a  i  X  A'  matrix,  whose  k-i\i  column  is  the  direction 
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vector  for  the  A;-th  signal,  and  r(/„)  =  [  ri(/„)     •  •  •     rA'(/„)  Y ^  is  the  K  x  I  signal 
vector,  and  N(/„)  =  [  NAfn)    ■  •  •     NL{fn)  1^'  ^^  the  I  x  1  noise  vector. 
The  resulting  spatial  covariance  matrix  is 

R(/„)  =  A(0,/„)R3(/„)A^(e,/„)  +  ^'(/n)R„(/n),  (2.84) 

where  R5(/„)  and  R„(/„)  are  the  signal  covariance  matrix  and  the  noise  covariance 
matrix,  respectively.  The  most  important  question  here  is  how  to  combine  the  in- 
formation from  different  frequency  bins.  The  idea  of  the  coherent  signal-subspace 
method  (CSM)  [29,  30]  is  to  use  transformation  matrices  to  focus  the  energy  of  the 
wideband  signals  to  a  single  reference  frequency.  The  focusing  removes  the  frequency 
dependency  of  the  data.  These  focusing  transformation  matrices  T(/n)  are  designed 
to  align  the  signal-subspace  at  all  frequencies  with  the  signal  subspace  at  the  reference 
frequency  /o, 

T(/„)A(0,/„)~A(0,/o),     n  =  !,•••,  J.  (2.85) 

Thus,  a  coherent  averaging  of  the  covariance  matrices  R(/„),  n  =  1,  ■  ■  • ,  J,  may  be 
performed  as  below 


R.       =       7  E  T(/n)R(/n)T^(/n) 
"'   n=l 

«    A(0,/o)R.A^(0,/o)  +  Rny  (2.86) 
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where 

1    -^ 
R,  =  -Y.^s{fn),  (2.87) 

■^   n=l 

and 

Rn.  =  7  E  a^(A)T(/„)R„(/„)T^(/„).  (2.88) 

^  n=l 

Hence  any  narrowband  subspace-based  algorithms  may  be  appHed  to  Rj,  after  the 
coherent  averaging  is  performed. 

It  is  apparent  from  (2.85)  that  the  transformation  matrix  depends  on  the  un- 
known AOAs  0.  To  form  the  focusing  transformation  matrices,  one  needs  to  have 
prehminary  estimates  of  AOAs.  The  preliminary  estimates  may  be  obtained  by  other 
methods,  such  as,  FFT  or  Capon's  method  [29,  30].  A  major  disadvantage  of  CSM 
is  that  the  focusing  may  cause  biased  estimates,  depending  on  the  accuracy  of  the 
preliminary  estimates.  We  will  proposed  a  wideband  RELAX  algorithm  for  angle 
estimation  of  wideband  plane  waves  arriving  at  a  uniform  linear  array.  The  wide- 
band RELAX  does  not  need  initial  estimates  as  other  wideband  algorithms  require 
and  naturally  focuses  the  narrowband  components  in  the  spatial  frequency  domain. 
Thus  the  wideband  RELAX  can  avoid  the  bias  problem  encountered  by  CSM  that 
uses  focusing  matrices. 

2.6     Polarization  Diversitv 

We  have  considered  an  array  of  omnidirectional  isotropic  sensors  so  far.  When 
a  polarization  sensitive  array  is  used,  the  output  of  sensor  array  is  related  to  the 
states  of  polarization  of  the  incident  plane  waves.  In  other  words,  ai{6k)  in  (2.26)  is 
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no  longer  constant,  it  depends  on  the  polarization  property  of  antenna  sensors  and 
incident  EM  plane  waves.  The  advantages  of  utilizing  polarization  sensitive  arrays 
have  been  discussed  previously  [10,  84,  85,  31,  86,  87,  33,  65,  88,  89,  90,  91].  In 
[10,  84]  the  MUSIC  algorithm  was  used  for  direction  finding  with  a  diversely  polarized 
antenna  array.  In  [92,  31,  86,  87]  the  methods  of  angle  and  polarization  estimation 
using  ESPRIT  with  a  crossed  dipole  array  were  presented.  A  noise  subspace  fitting 
method  with  diversely  polarized  antenna  arrays  was  proposed  in  [65].  The  discussion 
on  the  performance  analysis  of  diversely  polarized  antenna  array  can  be  found  in 
[88,  89].  A  ML  parameter  estimator  with  a  crossed  dipole  array  was  proposed  in 
[93]  for  partially  polarized  EM  waves.  However,  there  are  many  shortcomings  in  the 
previous  work.  For  example,  the  noise  subspace  fitting  and  MUSIC  cannot  work 
well  for  coherent  signals,  the  performance  of  the  crossed  dipole  array  is  sensitive  to 
AOAs.  New  approaches  will  be  proposed  in  Chapters  5  and  6  to  overcome  these 
problems.  We  will  consider  the  case  where  all  incident  narrowband  EM  plane  waves 
are  completely  polarized. 

2.7     Other  Applications 

The  research  progress  of  parameter  estimation  algorithms  in  sensor  array  pro- 
cessing has  resulted  in  a  great  diversity  of  applications,  such  as  multiuser  estimation 
and  detection,  channel  identification  and  spatial  diversity  in  personal  communica- 
tions. 

The  MUSIC  algorithm  was  used  in  [94]  for  propagation  delay  estimation  of  the 
DS-CDMA  (Direct-Sequence  Code  Division  Multiple  Access)  signals.    In  [95,  96]  a 
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large  sample  maximum  likelihood  (LSML)  algorithm  was  proposed  for  propagation 
delay,  carrier  phase,  amplitude  estimation  of  DS-CDMA  signals.  These  algorithms 
have  been  shown  to  be  robust  against  the  multiuser  near-far  problem  of  DS-CDMA 
systems.  Array  signal  processing  algorithms  are  expected  to  play  an  important  role 
in  accommodating  a  multiuser  communication  environment,  subject  to  severe  multi- 
path. 

2.8     Summary 

We  have  reviewed  important  properties  of  electromagnetic  waves  including  wave 
propagation  and  wave  polarization.  The  previous  work  in  sensor  array  processing 
has  been  reviewed  for  both  nonparametric  and  parametric  methods.  We  have  also 
described  the  relations  of  our  work  herein  with  the  previous  work.  A  rather  complete 
list  of  references  on  sensor  array  processing  can  be  found  in  [97]. 


CHAPTER  3 
NARROWBAND  ANGLE  AND  WAVEFORM  ESTIMATION  VIA  RELAX 


3.1     Introduction 

Many  high  resolution  array  processing  algorithms  have  been  devised  to  estimate 
the  incident  angles  of  the  signals  arriving  at  an  array  of  sensors.  The  success  of  those 
algorithms  is  typically  demonstrated  by  means  of  simulated  numerical  examples, 
often  under  the  additive  white  Gaussian  noise  assumption.  The  high  resolution  angle 
estimation  algorithms,  however,  are  usually  very  sensitive  to  the  violation  of  the  data 
model  assumptions.  As  a  result,  the  performance  of  those  algorithms  is  frequently 
much  worse  when  they  are  used  in  practical  applications.  To  derive  more  robust 
array  processing  algorithms,  we  shall,  for  instance,  relax  the  additive  white  noise 
assumption. 

In  this  chapter  we  describe  how  the  RELAX  algorithm,  recently  proposed  in  [45] 
for  mixed  spectrum  estimation,  can  be  used  for  angle  and  waveform  estimation  of 
narrowband  plane  waves  arriving  at  a  uniform  linear  array  in  the  presence  of  spatially 
colored  noise,  such  as  an  autogressive  (AR)  or  autogressive  moving  average  (ARMA) 
noise.  The  RELAX  algorithm  presented  in  this  paper  is  devised  for  the  multisnapshot 
case.  To  evaluate  the  performance  of  the  RELAX  algorithm,  we  apply  RELAX  to 
both  simulated  and  experimental  data,  and  compare  the  performance  of  RELAX 
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with  that  of  other  well-known  algorithms  including  ESPRIT  with  forward/backward 
spatial  smoothing  [12,  72],  MODE/WSF  [25,  23],  and  AP/ANPA  [66,  67]  under  the 
same  scenarios.  The  experimental  data  was  collected  by  the  Muti-parameter  Adaptive 
Radar  System  (MARS)  developed  at  McMaster  University.  More  importantly,  we 
explain  why  RELAX  outperforms  other  well-known  algorithms  including  ESPRIT, 
MODE/WSF,  and  AP/ANPA.  For  the  experimental  data,  since  we  do  not  know  the 
true  incident  angles,  we  introduce  a  cross-validatory  method  to  assess  the  quality  of 
the  estimates  obtained  with  the  algorithms  we  consider. 

3.2     Problem  Formulation 

Consider  the  problem  of  angle  and  waveform  estimation  of  K  narrowband  plane 
waves  impinging  on  a  uniform  linear  array  (ULA)  with  M  {M  >  2K)  elements. 
Assume  the  number  of  incident  plane  waves  K  is  already  known.  The  array  output 
vector  can  be  written  as: 

y(n)  =  A(0)s(n)  +  e(n),     n  =  1,2,  ■  ■  • ,  A^.  (3.1) 

In  (3.1),  N  is  the  number  of  temporal  snapshots,  e(n)  is  a  complex  M  x  1  noise 
vector,  s(n)  is  the  complex  K  x  1  signal  vector  ,  and 


A(0)  = 


a(0i)    a(^2)     •••    ^{Ok) 


(3.2) 


where  0 


01    02    ■■■    9 


K 


and  a(^fc),  k  =  1,2,  •  •  •  ,/^,  is  the  complex  M  x  1 


direction  vector  for  the  kih.  signal  arriving  from  angle  Ok  relative  to  the  array  normal. 
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(Here  (•)-^  denotes  the  transpose.)   For  the  ULA,  the  kth  direction  vector  a(^yt)  has 
the  form: 


^{e,)  = 


-iT 


1     e^'^**'"^^ 


gi2fi(M-l)sin0fc 


,      k  =  l,2,---,K, 


(3.3) 


where  6  is  the  spacing  between  the  array  sensors  and  A  is  the  signal  wavelength. 

The  noise  vector  e(n)  is  a  spatially  correlated  random  process  that  is  assumed 
to  be  temporally  white.  In  particular,  e(n)  has  zero-mean  and  unknown  covariance 


matrix  Q: 


E{e(ni)e^(n2)}  =  QSn,,n2, 


(3.4) 


where  (•)^  denotes  the  complex  conjugate  transpose  and  <5„j,„2  i^  ^^e  Kronecker  delta. 

T 


The  signal  waveforms  s(n) 


,  n  =  1,2,---,A^, 


•si(n)    S2{n)    ■■■    SK{n) 
are  modeled  as  deterministic  unknowns.    This  assumption  is  usually  referred  to  as 
the  deterministic  (or  conditional)  signal  model  [20,  76]. 

The  problem  of  interest  herein  is  to  estimate  the  angles  6i,  92,  •  •  •,  Ok  and  the 
complex  waveforms  5i(n),  S2{n),  •  •  •,  ^^-(n)  from  y(n),  n  =  1,2,  •  •  • ,  A^. 

3.3     Angle  and  Waveform  Estimation  Using  RELAX 

The  RELAX  algorithm  was  derived  in  [45]  as  an  estimator  of  sinusoidal  param- 
eters in  the  presence  of  colored  noise.  We  describe  below  how  the  RELAX  algorithm 
can  be  modified  for  angle  and  waveform  estimation  with  a  uniform  linear  array  of 
sensors. 
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The  angle  and  waveform  estimates  {^,  s(l),  s(2),  ■  •  • ,  s(A'')}  can  be  obtained  via 
RELAX  by  minimizing,  with  respect  to  0  and  s(l),  s(2),  •••,  s(A^),  the  following 
nonlinear  least-squares  (NLS)  criterion: 


F,  [0,  s(l),  s(2),  •  •  • ,  s{N)]  =  f^Mn)  -  A{0)s{n)]"[y{n)  -  A{0)s{n)].         (3.5) 

n=l 


Note  that  if  the  noise  e(n)  were  spatially  white  and  Gaussian,  minimizing  (3.5)  would 
have  yielded  the  deterministic  maximum  likelihood  estimates  of  the  incident  angles 
and  signal  waveforms  [28,  66]. 

We  present  below  our  approach  to  the  minimization  of  the  cost  function  in  (3.5). 
Assume  for  now  that  there  are  K  signals,  where  K  is  an  intermediate  value  (  K  <  K, 
see  the  steps  below).  Let 


yk{n)^y{n)-     J2    a(^0-5t(n),  (3.6) 


where  {di,Si{l),Si{2),-  ■  ■  ,Si{N)}^^^-^f,  are  assumed  to  have  been  previously  esti- 
mated. Then  the  cost  function  for  estimating  the  kth  signal  parameters  becomes 

N 

F2{0k,  6,(1),  .,(2),  •  •  • ,  Sk{N))  =  Y:  [ykin)  -  Bi{ek)s,{n)f  [y,(n)  -  a(^,).fc(n)] , 

n=l 

(3.7) 
The  minimization  of  F2  with  respect  to  6k  and  Sk{l),-  ■  ■ ,  Sk{N)  gives 


.  a^(^fc)yfc(n) 

Sk{n)  =  - 


M 


_  ,     n  =  l,2,---,N,  (3.8) 
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and 


N 


0k  =  arg  min  ^ 


n=l 


I- 


a(^;c)a^(^fc) 


M 


yk{n) 


N 


argmax^  a^(^fc)y;t(n)     .      (3.9) 


n-l 


Hence  6k  is  obtained  as  the  location  of  the  dominant  peak  of  the  sum  of  the  peri- 

2 

odograms  a^(^fc)yfc(n)  /M,  over  n  =  1, 2,  •  •  • ,  A^,  which  can  be  efficiently  computed 
by  using  FFT  (fast  Fourier  transform)  with  each  of  the  data  sequence  Ykin)  padded 
with  zeros.  (Note  that  padding  with  zeros  is  necessary  to  determine  6k  with  high 
accuracy.)  Furthermore,  Sk{n)  is  easily  computed  from  the  complex  height  of  the 
peak  of  a."{6k)yk{n)/M  (for  n  =  1, 2,  •  •  • ,  iV). 

We  can  now  proceed  to  present  our  relaxation  (RELAX)  algorithm  for  the  min- 
imization of  the  nonlinear  least-squares  cost  function  in  (3.5)  .  Note  that  we  assume 
the  number  of  incident  signals  K  is  already  known  or  estimated.  The  RELAX  algo- 
rithm comprises  the  following  steps: 


Step   (1):    Let  K  =   1.    Obtain  6i  and  51(77)  from  yi(n)  as  in  (3.8)  and  (3.9), 
n  =  l,2,---,Af. 


Step  (2):  Let  K  =  2.  Compute  y2(^)  with  (3.6)  by  using  9i  and  Si{n)  obtained  in 
Step  (1),  n  =  1,2,-  ■  ■  ,N.  Obtain  62  and  S2(n)  from  y2{n)  as  in  (3.8)  and  (3.9), 
n  =  1,2,---,A^. 

Next,  compute  yi(n)  with  (3.6)  by  using  62  and  S2{n)  and  redetermine  61  and 
si{n)  from  yi(n)  as  in  (3.8)  and  (3.9),  n  =  1,2,-  ■  •  ,N. 
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Iterate  the  previous  two  substeps  until  "practical  convergence"  is  achieved  (to 
be  discussed  later  on). 

Step  (3):  Let  K  =  3.  Compute  y3(n)  with  (3.6)  by  using  {4-,5,-(l),5i(2),  ■  •  ■  ,^N)}f^i 
obtained  in  Step  (2).    Obtain  63  and  S3(n)  from  ysin)  as  in  (3.8)  and  (3.9), 
n  =  1,2,---,N. 

Next,  compute  yi(n)  with  (3.6)  by  using  {^i,  s,(l),  5,(2),  •  •  ■  ,Si{N)}^^2  ^-nd  re- 
determine §1  and  si{n)  from  yi(n),  n  =  1,2,-  •  •  ,N.  Then  compute  y2(n)  with 
(3.6)  by  using 

{0i,Si{l),Si{2),-  ■  ■  ,Si{N)}i=i^3  and  redetermine  62  and  S2{n)  from  y2(n),  n  = 
1,2,---,7V. 

Iterate  the  previous  three  substeps  until  "practical  convergence"  is  achieved. 

Remaining  Steps:  Continue  similarly  until  K  is  equal  to  the  known  number  K  or 
the  estimated  number  K  of  signal  sources.  (See  [45]  for  a  possible  approach  to 
the  estimation  of  K  when  it  is  unknown.) 

The  "practical  convergence"  in  the  iterated  substeps  of  the  above  RELAX  algo- 
rithm may  be  determined  by  checking  the  relative  change  of  the  cost  function  in  (3.5) 
between  two  consecutive  iterations.  In  our  simulated  and  experimental  examples,  to 
be  presented  in  the  next  subsection,  we  terminate  the  iterations  when  the  relative 
change  is  less  than  or  equal  to  e  =  10~^,  i.e.. 


|Fi[g,+i,s,4.i(l),---,s,+i(JV)]-Fi[g„s,(l),---,s,(iV)]|  ^^  ^^^^^ 


F,[e„s,{l),---,s,{N)] 
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where  ^,,8,(1),  •  ■  •  ,Sq{N)  denote  the  estimates  obtained  after  the  9th  iteration. 

To  show  how  RELAX  works,  we  consider  an  example  shown  in  Figure  3.1,  where 
A^  =  1,  K  =  2,61  =  35.5°,  62  =  38.5°,  si{l)  =  1,  62(1)  =  1,  and  the  ULA  has  M  =  32 
sensors  with  the  spacing  between  two  adjacent  sensors  equal  to  a  half  wavelength  of 
the  incident  plane  waves.  The  frequency  axis  used  in  the  figure  is  related  to  AOAs 
by  fk  =  |sin(^fc),  ^  =  1,2,  where  /i  and  /2  are  shown  with  the  dashed  Hue.  The 
FFTs  of  Yk  in  Figure  3.1  is  obtained  when  the  array  data  are  zero  padded  to  1024. 
Note  that  the  AOA  separation  3°  is  less  than  the  FFT  resolution  3.6°.  However,  the 
AOAs  can  be  estimated  accurately  with  RELAX  after  a  few  of  iterations. 

3.4     Results  of  Simulated  and  Experimental  Data 

We  first  present  a  simulation  example  for  estimating  the  angles  and  waveforms 
of  incident  signals.  The  simulation  results  were  obtained  by  using  50  Monte-Carlo 
trials.  The  noise  em(n)  is  assumed  to  be  a  spatial  AR  random  process  of  order  1,  i.e., 

em{n)^ -aiem-i{n)  +  Wm{n),  (3.11) 

where  Oi  =  —0.85  and  Wm{n),  m  =  1,2,  •••,M,  are  independently  and  identically 
distributed  zero-mean  complex  Gaussian  random  variables  with  variance  cr^.  The 
complex  amplitudes  of  the  waveforms  are  assumed  to  be  1.  The  SNR  referred  to 
below  is  defined  by  [98] 


SNR=101ogiJ^J^j    [dB].  (3.12) 
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We  assume  that  there  are  K  —  2  correlated  incident  signals,  with  the  correlation 
coefficient  equal  to  0.99,  arriving  from  angles  ^i  =  —10°  and  02  =  2°.  The  number  of 
incident  signals  is  assumed  known.  The  array  is  assumed  to  have  M  =  8  sensors  that 
axe  uniformly  spaced  with  the  spacing  between  two  adjacent  sensors  equal  to  a  half 
wavelength.  The  number  of  snapshots  taken  at  the  array  output  is  A'^  =  10.  All  array 
output  vectors  are  zero  padded  to  2^^  before  using  them  with  FFT  in  the  RELAX  and 
ANPA  algorithms.  (Note  that  the  higher  the  SNR,  the  more  padded  zeros  are  needed 
so  that  the  estimation  accuracy  is  not  affected  by  the  FFT  intervals.)  The  waveform 
estimates  for  MODE,  ANPA,  and  LS-ESPRIT  (least  squares  based  ESPRIT)  [12]  are 
obtained  by  using 

s(n)  =  [A"{e)A{e)]-'A"{e)y{n).  (3.13) 

The  root-mean-squared  error  (RMSE)  of  the  ^th  waveform  estimate  is  defined  as 


y^n=i  \^k{n)  —  ifc(n)p.  Note  that  LS-ESPRIT  is  used  with  forward/backward  spa- 
tial smoothing  and  subarray  length  equal  to  4  (that  is,  M/2)  due  to  the  highly 
correlated  signals. 

Figures  3.2(a)  and  (b)  show  the  RMSEs  of  the  angle  and  waveform  estimates, 
respectively,  for  the  first  signal  as  a  function  of  SNR.  (The  results  for  the  second  signal 
are  similar.)  The  estimator  performances  are  also  compared  to  the  corresponding 
Cramer-Rao  bound  (CRB).  (The  CRB  matrix  for  the  signal  parameters  is  derived  in 
Appendix  A.)  Note  that  RELAX  has  the  best  performance  among  the  four  algorithms 
for  both  angle  and  waveform  estimation,  especially  for  waveform  estimation  at  low 
SNR.  RELAX  outperforms  MODE  and  ESPRIT  since  the  noise  is  not  spatially  white 
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as  the  latter  methods  assume.  ANPA  gives  bad  waveform  estimates  at  low  SNR 
because  using  ANPA  may  yield  almost  identical  angle  estimates  in  some  Monte- 
Carlo  trials.  The  ANPA  angle  estimates  are  obtained  by  maximizing  the  following 
cost  function: 

^3W  =  EIIPA(0)yHf'  (3-14) 


n=l 


which  can  be  obtained  by  first  minimizing  the  cost  function  in  (3.5)  with  respect  to 
s(l),s(2),---,s(A^)  [66,  67],  where  Pa  =  a(a^a)~^  A^.  Although  using  (3.13) 
for  waveform  estimation  assumes  that  the  angle  estimates  are  not  nearly  identical, 
maximizing  F3{0)  does  not  guarantee  this  assumption.  Figure  3.3  shows  F3{G)  as  a 
function  of  9i  and  02  when  the  realization  of  the  noise  is  the  one  that  gives  the  worst 
waveform  estimates  for  ANPA  in  Figure  1.  It  can  be  seen  from  Figure  3.3  that  the 
maximum  of  the  cost  function  corresponds  to  almost  identical  angles.  Hence  maxi- 
mizing the  cost  function  in  (3.14)  with  respect  to  6  (only)  and  using  the  estimated 
0  to  estimate  s(n)  with  (3.13),  instead  of  leading  to  a  simpler  problem,  can  actually 
complicate  the  optimization  of  the  cost  function  in  (3.5)  and  lead  to  poorer  estimates 
when  the  incident  angles  are  closely  spaced. 

We  next  apply  the  algorithms  to  the  experimental  data  collected  by  the  array 
system  known  as  the  Multi-parameter  Adaptive  Radar  System  (MARS)  [99].  The 
array  system  was  developed  at  the  Communications  Research  Laboratory  at  McMas- 
ter  University.  The  data  was  collected  by  deploying  the  array  at  the  west  coast  of 
the  Bruce  Peninsula,  Ontario,  Canada,  overlooking  Lake  Huron.  MARS  is  a  vertical 
uniform  linear  array  consisting  of  M  =  32  horizontally  polarized  horn  antennas.  The 


50 


spacing  between  adjacent  antenna  sensors  is  5.715  cm.  The  four  sets  of  data  we  use 
below  were  collected  when  the  array  system  was  operated  at  frequencies  9.76,  9.79, 
11.32  and  12.34  GHz.  The  data  was  recorded  with  r2-bit  precision  and  sampled  at 
62.5  samples  per  second.  For  each  carrier  frequency,  127  snapshots  were  collected  at 
each  antenna  output.  There  are  two  incident  signals,  a  fact  assumed  to  be  known  to 
the  angle  estimation  algorithms.  One  of  the  incident  signals  is  the  direct  path  and 
the  other  is  the  specular  path,  which  is  reflected  from  the  lake.  The  direct  incident 
signal  is  a  continuous  wave  (CW),  whose  amplitude  is  a  constant  and  whose  phase  is 
a  linear  function  of  time.  Since  the  specular  path  is  a  delayed  and  reflected  version 
of  the  direct  path,  the  phase  difference  between  the  two  paths  is  a  constant,  which 
is  determined  by  the  time  delay  between  the  direct  and  specular  paths,  the  carrier 
frequency  of  the  waves,  and  the  reflection  coefficient  of  the  lake.  The  incident  signals 
arrive  from  near  the  array  normal,  but  the  exact  incident  angles  are  unknown  since 
the  vertical  array  structure  may  have  been  on  a  slight  tilt.  The  parameter  estimation 
algorithms  we  consider  below  do  not  assume  any  a  priori  knowledge  of  the  incident 
angles  and  the  signal  waveforms. 

In  this  application  we  use  a  cross-validatory  criterion  to  assess  the  quality  of 
estimates.  Consider  the  case  of  where  the  algorithms  are  applied  to  a  single  snapshot 
at  a  time.  Let  9s  be  the  estimated  angle  vector  obtained  from  the  5th  snapshot  by 
using  one  of  the  algorithms.  Then  our  cross-validatory  criterion  is: 


V  /   s=l  n=l,n:^5 
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where  N  is  the  number  of  those  snapshots  that  do  not  give  almost  identical  angle 
estimates,  and  P^  =  I— A  (A^AJ  A^.  (Note  that  almost  identical  angle  estimates 
cause  the  matrix  A^A  in  P-"-  ;,  to  be  ill  conditioned.)  What  the  criterion  does  here 
may  be  explained  by  the  following  steps:  1)  estimate  AOAs  from  a  single  snapshot;  2) 
obtain  least  squares  estimate  of  waveforms  by  using  the  estimate  of  AOAs  for  every 
other  snapshot,  i.e.,  s(n)  =  [A^(^,)A(^)]-iA^(^s)y(n),(n  =  1,  •  •  •  ,iV,  and  rz  7^  s); 
3)  obtain  the  residue  by  subtracting  the  signal  part  from  the  data  with  the  data  model 
for  each  snapshot,  i.e,  residuehi)  =  y(n)  —  A{Os)s{n)  =  P-*-  -  y(n),(ra  =  1,  •  •  • ,  A'', 
and  n  /  5);  4)  find  the  summation  of  residues  for  each  of  these  snapshots  as  in  (3.15). 
Clearly,  the  smaller  the  summation,  the  better  the  performance  of  the  algorithm.  The 
inner  sum  in  (3.15)  shows  how  well  (or  bad)  the  estimate  obtained  from  snapshot  s 
can  be  used  to  predict  the  array  outputs  observed  in  the  other  available  snapshots, 
according  to  the  assumed  data  model.  The  outer  sum  adds  those  prediction  errors 
together.  The  lower  the  C,  the  better  the  performance  of  the  algorithm. 

Table  1  shows  the  cross- validatory  criterion  C  obtained  by  using  RELAX,  ANPA, 
MODE,  LS-ESPRIT,  and  TLS-ESPRIT  (total  least  squares  based  ESPRIT)  [12]  at 
frequencies  8.62,  9.76,  9.79,  and  12.34  GHz.  Note  that  LS-ESPRIT  and  TLS-ESPRIT 
are  used  with  forward/backward  spatial  smoothing  and  subarray  length  equal  to  10. 
In  Table  1,  the  failure  rate  indicates  how  often  an  algorithm  gives  almost  identical 
angle  estimates,  and  hence  very  poor  waveform  estimates.  Note  that  for  the  data  set 
analyzed  here,  RELAX  and  ANPA  have  similar  performances  and  are  better  than 
the  eigendecomposition  based  MODE  and  ESPRIT,  which  assume  that  the  additive 
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noise  is  spatially  white  (note  that  the  noise  statistics  are  unknown).  To  see  whether 
the  noise  is  actually  spatially  correlated,  we  estimate  the  noise  covariance  matrix  as 
follows: 

Figure  3.4  shows  the  correlation  coefficients  between  the  first  sensor  noise  and  the 
other  sensor  noises.  The  noise  appears  to  be  strongly  spatially  correlated,  which 
explains  the  poorer  performance  of  the  eigenstructure  based  methods.  Finally,  we 
note  from  Table  1  that  LS-ESPRIT  and  TLS-ESPRIT  have  similar  performances  for 
this  data  set,  but  TLS-ESPRIT  has  a  much  higher  failure  rate,  which  indicates  that 
when  waveform  estimation  is  desired,  LS-ESPRIT  should  be  preferred  over  TLS- 
ESPRIT. 

Finally,  we  compare  the  computational  complexities  of  these  algorithms.  The 
amount  of  computations  needed  by  RELAX  depends  on  N,  K,  SNR  and  the  num- 
ber of  zero-paddings.  For  example,  RELAX  requires  about  40%  of  the  amount  of 
computations  required  by  ANPA  when  N  =  I,  K  =  2,  SNR=  0  dB  and  the  num- 
ber of  zero-paddings  is  2^^.  Since  ESPRIT  and  MODE  do  not  need  iterations,  they 
require  much  fewer  computations  than  RELAX  and  ANPA.  For  example,  ESPRIT 
and  MODE  only  need  about  0.3%  and  1.5%,  respectively,  of  the  amount  of  compu- 
tations required  by  RELAX  when  iV  =  1,  A'  =  2,  SNR=  0  dB  and  the  number  of 
zero-paddings  is  2^^.  However,  RELAX  can  still  be  attractive  because  it  can  be  easily 
implemented  with  simple  EFT  chips  in  parallel. 
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3.5     Conclusions 

We  have  presented  a  RELAX  algorithm  for  the  multi-snapshot  case  of  angle  and 
waveform  estimation  of  narrowband  plane  waves  arriving  at  a  uniform  linear  array. 
The  RELAX  algorithms  are  conceptually  and  computationally  simple;  their  imple- 
mentations mainly  require  a  sequence  of  fast  Fourier  transforms.  We  have  evaluated 
the  performance  of  the  RELAX  algorithm  by  using  both  simulated  and  experimental 
data,  and  compared  it  with  the  performances  of  other  well-known  algorithms  such  as 
ESPRIT  with  forward/backward  spatial  smoothing,  MODE/WSF,  and  AP/ANPA. 
We  have  explained  by  means  of  results  of  both  simulated  and  experimental  data  why 
better  signal  parameter  estimates  can  be  obtained  by  using  RELAX  as  compared  to 
using  the  other  algorithms. 


54 


Table  1:  Cross- validatory  criterion  for  the  different  algorithms  when  used  with  the 

experimental  data  collected  by  MARS. 


Algorithms 

8.62  GHz 

9.76  GHz 

9.79  GHz 

12.34  GHz 

RELAX 

Criterion  C 

1.3578e+04 

2.9153e+04 

6.6079e+04 

1.7177e+04 

RELAX 

Failure  rate 

0 

0 

0 

0 

ANPA 

Criterion  C 

1.3323e+04 

2.9209e+04 

6.6260e+04 

1.7172e+04 

ANPA 

Failure  rate 

0 

0 

0 

0 

MODE 

Criterion  C 

4.5306e+04 

3.8840e+04 

7.2059e+04 

2.6091e+04 

MODE 

Failure  rate 

38% 

9% 

5% 

6% 

LS- 
ESPRIT 

Criterion  C 

1.5353e+04 

3.3716e+04 

11.0460e+04 

1.8302e+04 

LS- 
ESPRIT 

Failure  rate 

0 

0 

0 

0 

TLS- 
ESPRIT 

Criterion  C 

1.6723e+04 

3.3712e+04 

11.4600e+04 

1.8302e+04 

TLS- 
ESPRIT 

Failure  rate 

47% 

0 

9% 

0 
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Figure  3.1.  An  example  of  using  RELAX,  (a)  Modulus  of  FFT  of  y  in  step  (1), 
K  =  1.  (b)  Modulus  of  FFT  of  y2  in  step  (2),  /?  =  2  (1st  iteration),  (c)  Modulus  of 
FFT  of  yi  in  step  (2),  K  =  2  (2nd  iteration),  (d)  Modulus  of  FFT  of  ya  in  step  (2), 
K  =  2  (3rd  iteration),  (e)  Modulus  of  FFT  of  yi  in  step  (2),  /?  =  2  (4th  iteration), 
(f)  Modulus  of  FFT  of  y2  in  step  (2),  K  =  2  (5th  iteration). 
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(b) 

Figure  3.2.  RMSEs  of  the  angle  and  waveform  estimates  of  the  first  signal  as  a 
function  of  SNR  when  iV  =  10,  M  =  8,  and  K  =  2  correlated  incident  signals  with 
the  correlation  coefficient  equal  to  0.99  arrive  from  9i  =  -10°  and  O2  =  2°.  (a)  Angle 
estimation,  (b)  Waveform  estinaation. 
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(b) 

Figure  3.3.  F3  vs.  9i  and  62  when  A'^  =  10,  M  =  8,  the  correlation  coefficient  of 
the  two  incident  signals  is  equal  to  0.99,  and  the  realization  of  the  noise  is  the  one 
that  gives  the  worst  waveform  estimates  for  ANPA  in  Figure  1.  (a)  Mesh  plot,  (b) 
Contour  plot. 
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Figure  3.4.  Estimated  noise  correlation  coefficients  between  the  first  and  the  other 
sensors.  Figures  (a)  -  (d)  are  for  the  carrier  frequencies  8.62,  9.76,  9.79,  and  12.34 
GHz,  respectively. 


CHAPTER  4 
ANGLE  ESTIMATION  OF  WIDEBAND  SIGNALS  USING  RELAX 

4.1     Introduction 

In  this  chapter  we  extend  the  RELAX  algorithm  to  the  case  of  wideband  sources. 
We  concentrate  on  angle  estimation  herein.  We  show  that  the  wideband  RELAX 
algorithm  we  devise  naturally  focuses  the  narrowband  components  in  the  spatial 
frequency  domain.  We  use  both  numerical  and  experimental  examples  to  demonstrate 
the  performance  of  the  wideband  RELAX  algorithm  and  compare  its  performance 
with  that  of  the  well-known  CSM  [29]  based  LS-ESPRIT  (CSM-ESPRIT).  We  also 
explain  why  wideband  RELAX  can  outperform  the  CSM-ESPRIT. 

4.2     Problem  Formulation 

Consider  the  same  uniform  Hnear  array  (ULA)  with  M  elements  as  before.  As- 
sume that  the  incident  wideband  deterministic  signals  have  a  common  bandwidth 
B  (Hz)  with  center  frequency  /o  (Hz).  The  kth.  bandpass  signal  Sk{t)  observed  at  a 
reference  point  can  be  written  as 

Skit)  =  7fc(0e^'^^°',  (4-1) 
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where  "ykit)  denotes  the  complex  envelope.  Let  the  signal  be  observed  over  a  duration 
[to,to  +  To].  The  complex  envelope  can  be  written  as 


7k{t)  =  E  rfc(/;)e''^^'\     to<t<to  +  To,  (4.2) 


where 

1      fto+To  .^    .  ^        ^ 

Ufi)  =  ^  ^k{t)e-^'-^^'dt,  (4.3) 

IQ  J  to 

with  //  =  '~<^+^)/^^  /  =  1,2,  •••,i.  L  is  the  number  of  frequency  components 
symmetrically  placed  around  0  Hz  with  /;+i  -  /;  =  j^  =  ^.  Thus  the  bandpass 
signal  at  the  reference  point  (say  the  first  sensor)  can  be  written  as 


L 

I 

1=1 


Sk{t)  =  J2r,ifi)e^'^^^°^^'^'.  (4.4) 


The  jfcth  signal  at  the  mth  sensor  has  a  propagation  time  delay  n%,„  so  that 


L 

I 

1=1 


s,{t  +  n,m)  =  ^r,(/Oe^'^-(^''+^')('+^^-)  =  x,{t,n,^)e^'^^°\  (4.5) 


where 

(m-l)6   . 
Tk,m  = sm 


Ok,     k  =  l,2,---,K,     m  =  l,2,---,M,  (4.6) 


C 

with  8  being  the  spacing  between  sensors,  C  the  propagation  speed,  and  ^fc  the  angle 
of  arrival  of  the  A;th  incident  signal  relative  to  the  array  normal. 
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Let  Xfc(i)  = 


Xk{t,Tk,l)      Xk{t,Tk,2)      ■•■      Xk{t,TkM) 


.  Then 


1=1 


j2,vf,t 


(4.7) 


where  a{9k,  /o  +  //)  is  the  M  x  1  direction  vector  of  the  fcth  source  and  has  the  form 


a{9kJo  +  fi) 


-iT 


oJ27r(/o+/i)Tfc,i       pj27r(/o+/()rfc,2       .  .  .      pj'2-^(fo+fl)rk,M 


(4.8) 


Hence  if  K  wideband  signals  along  with  some  noise  simultaneously  impinge  upon 
the  sensor  array,  the  received  data  vector  has  the  form 


L 

I 


x(0  =  Y:  [A(e,  /o  +  //)r(/)  +  N(/)]  e^'-^"  =  X:  y(/)e^''^^'S  (4.9) 


/=i 


where 


M^Jo  +  fi) 


a(^i,/o  +  //)    •■■    a(^A-,/o  +  /0 


(4.10) 


and 


r(/) 


ri(/0  •••  TKifi) 


is  the  K  X  1  signal  vector,  and 


(4.11) 


N(0 


A^i(//)    ••■    NmUi) 


(4.12) 
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is  the  M  X  1  noise  Fourier  coefficient  vector,  and  {y(/)}  are  by  definition  the  M  x  1 
Fourier  coefficient  vectors  of  x(t): 

y(/)  =  A(0,/o  +  //)r(/)  +  N(/),     /=1,2,---,L.  (4.13) 

We  assume  that  the  noise  vector  N(/)  has  zero-mean  and 

E{N(/i)N^(/2)}  =  Q<^/„,„  (4.14) 

where  Q  is  unknown. 

The  problem  of  interest  herein  is  to  estimate  ^i,  62-,  ■■■■,&K  from  y(/),  /  = 
1,2,  •••,X.  The  main  difference  with  respect  to  the  narrowband  problem  treated 
in  Section  2  is  that  now  the  array  transfer  matrix  A  depends  on  the  snapshot  index 
/. 

4.3     Angle  Estimation  Using  RELAX 

The  estimates  {0, r(l),  •  •  •  ,T{L)}  can  be  obtained  via  RELAX  by  minimizing 
the  following  nonlinear  least-square  criterion  with  respect  to  0  and  r(l),---,  T{L) 
(See  (3.5)): 


Gi  =  E  [y(0  -  A(^,  /o  +  //))r(0]''  [y(/)  -  A(^,  /o  +  fi))r{i)] .         (4.i5) 
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To  estimate  the  parameters  of  the  A;th  signal,  consider 


1=1 


G2  =  E  MO  -  ^i^k,  /o  +  //))r.(0]''  [yk{i)  -  a(^,,  /o  +  //))r,(/)] ,       (4.16) 


where 

yfc(0  =  y(0-  E  ^{kfo  +  miil  (4.17) 

and  where  K  is  as  defined  before,  and  {6i,  Ti{l),  •  •  • ,  ^i{L)}^^  -^f.  are  assumed  to  be 
given  or  estimated.  Let 


Ml)  =  ^^(^^-/°^/')y.(0^  (4_^g) 


which  can  be  obtained  by  using  FFTs  with  zero  padding.   Then  similar  to  the  nar- 
rowband RELAX  algorithm,  Ok  can  be  obtained  by 


L 

^fc  =  argmax^|n(Ol',  (4.19) 


and 

f,(/)=n(/)|,^^^^,     l  =  l,2,...,L.  (4.20) 

Hence  6k  is  obtained  from  the  location  of  the  dominant  peak  of  the  sum  of  the 
focused  periodograms  |Fa:(/)|^  ,  /  =  1, 2,  •  •  • ,  L.  Then  tk{l)  is  easily  computed  from 
the  complex  height  of  the  peak  of  1^(0-  Note  that  the  focusing  of  the  narrowband 
components  can  be  naturally  achieved  in  (4.19)  by  expanding  and  compressing  the 
FFTs  (with  zero  padding)  of  y;t(/)  according  to  the  ratio  (/o  +  //)//o  since  the  spacing 
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between  two  adjacent  FFT  samples  are  different  for  different  /  .    The  steps  of  our 
RELAX  algoritfim  for  minimizing  (6.20)  are  summarized  as  follows: 

Step  (0):  Obtain  y(/)  from  x(^)  via  DFT  (discrete  Fourier  transform). 

Step  (1):  Assume  K  =  1.  Obtain  ^i  and  fi(/)  from  y(/)  as  in  (4.19)  and  (4.20), 
/=  1,2,---,L. 

Step  (2):  Assume/?  =  2.  Compute y2(/)  with  (4.17)  by  using  ^i  and  ri(/)  obtained 
in  Step  (1),  /  =  1,2,  ■••,!.  Obtain  §2  and  f2(0  from  y2(/)  as  in  (4.19)  and 
(4.20), /  =  1,2,  •••,!. 

Next,  compute  yi(/)  with  (4.17)  by  using  02  and  r2(/)  and  redetermine  61  and 
fi(Ofromyi(/), /=1,2,---,L. 

Iterate  the  previous  two  substeps  until  "practical  convergence"  is  achieved.  (See 
Section  2  for  details.) 

Remaining  Steps:  Continue  similarly  until  K  is  equal  to  the  known  number  K  or 
the  estimated  number  K  of  signals. 

We  refer  to  this  algorithm  as  the  wideband  RELAX. 

4.4     Numerical  Results 

We  present  below  numerical  and  experimental  examples  showing  the  perfor- 
mance of  the  proposed  algorithm  for  estimating  the  incident  angles  of  wideband 
signals  in  the  presence  of  either  white  or  unknown  colored  noise.  We  also  compare 
the  performance  of  wideband  RELAX  with  that  of  CSM  based  LS-ESPRIT  (CSM- 
ESPRIT)  for  angle  estimation. 
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In  both  of  the  simulation  examples  below,  the  array  is  assumed  to  be  a  ULA  of 
M  =  8  sensors  with  the  spacing  between  two  adjacent  sensors  equal  to  half  of  the 
wavelength  corresponding  to  the  center  frequency  /q.  The  wideband  sources  have 
the  same  center  frequency  /o  =  100  Hz  and  the  same  bandwidth  B  =  AO  Hz.  The 
noise  is  a  temporally  stationary  zero-mean  white  Gaussian  process,  independent  of 
the  signals,  and  spatially  either  a  white  or  an  autoregressive  (AR)  process.  The  total 
observation  time  is  To  =  0.8  seconds.  The  demodulated  data  is  sampled  at  twice  the 
Nyquist  rate.  The  array  output  is  decomposed  into  L  =  33  narrowband  components 
via  the  DFT  (discrete  Fourier  transform).  The  signal-to-noise  ratio  (SNR)  is  defined 
as  the  ratio  of  the  power  of  each  signal  to  the  noise  power.  We  assume  that  there 
are  A'  =  2  independent  wideband  signals  with  equal  power  impinging  on  the  array 
from  6i  =  10°  and  62  =  20°.  Note  that  the  two  signals  cannot  be  resolved  by 
the  spatial  periodogram.  All  sequences  are  zero  padded  to  4096  for  FFT  in  the 
RELAX  algorithm.  The  simulation  results  were  obtained  by  using  30  Monte-Carlo 
simulations. 

Figure  4.1(a)  shows  the  root-mean-squared  error  (RMSE)  of  the  second  signal 
as  a  function  of  SNR  in  the  presence  of  white  noise.  (The  results  for  the  first  signal 
are  similar.)  We  note  that  the  RMSE  of  the  wideband  RELAX  decreases  with  SNR, 
but  the  RMSE  of  CSM-ESPRIT  stops  decreasing  when  SNR  reaches  a  certain  level. 
This  is  because,  due  to  focusing,  CSM-ESPRIT  provides  biased  angle  estimates  even 
when  the  angle  separation  of  the  two  sources  is  within  the  resolution  of  the  spatial 
periodogram. 
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In  Figure  4.1(b),  we  consider  the  performance  of  both  wideband  RELAX  and 
CSM-ESPRIT  in  the  presence  of  unknown  AR  noise.  The  noise  Nm{fi)  is  assumed 
to  be  a  complex  AR  process  of  order  1,  i.e,  Nm{fi)  =  — ai-/Vm-i(//)  +  ^mift),  where 
Gi  =  — 0.85e~-'^  and  Wm{fi)  is  a  zero-mean  complex  white  Gaussian  random  process 
with  variance  equal  to  cr^.  The  SNR  is  defined  by  —  lOlog^Q  ^T  p  dB.  Note  that 
wideband  RELAX  again  performs  better  than  CSM-ESPRIT. 

Finally,  we  apply  both  RELAX  and  CSM-ESPRIT  to  the  experimental  data 
collected  by  the  sensor  array  testbed  [100,  101]  at  the  University  of  Minnesota.  The 
uniform  linear  array  consists  of  8  sensors.  The  spacing  between  adjacent  sensors  is 
about  2.1  times  the  wavelength  corresponding  to  /q.  Two  correlated  sources  arrive 
from  around  —33°  and  —36°.  The  SNR  is  21  dB  for  each  source.  The  center  carrier 
frequency  for  this  data  is  40  kHz.  The  bandwidth  of  the  data  is  4  kHz.  The  data  was 
sampled  at  a  rate  of  5  kHz,  and  was  decomposed  into  L  =  5  frequency  bins  (38.125 
kHz,  39.0625  kHz,  40  kHz,  40.9375  kHz,  and  41.875  kHz).  Figure  4.2  shows  the  angle 
estimates  obtained  from  several  observation  intervals.  (Note  that  one  of  the  angle 
estimates  does  not  show  up  in  Figure  4.2  (b)  because  it  is  too  small.)  The  means 
and  the  standard  deviations  of  the  angle  estimates  in  Figure  4.2  were  calculated 
by  averaging  the  angle  estimates  obtained  from  all  64  observation  intervals.  Note 
that  wideband  RELAX  provides  smaller  standard  deviations  and  biases  than  CSM- 
ESPRIT. 
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4.5     Conclusions 

A  wideband  RELAX  algorithm  for  the  angle  estimation  of  wideband  sources 
has  been  presented.  The  wideband  RELAX  naturally  focuses  the  narrowband  com- 
ponents in  the  spatial  frequency  domain.  Numerical  and  experimental  examples  have 
shown  that  the  wideband  RELAX  can  perform  better  than  CSM. 
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(a) 


(b) 


Figure  4.1.  RMSEs  of  the  angle  estimates  of  the  second  signal  as  a  function  of  SNR 
when  K  =  2  uncorrelated  wideband  signals  arrive  from  9i  =  10°  and  O2  =  20°, 
M  =  8,  and  L  =  33.  (a)  In  the  presence  of  white  noise,  (b)  In  the  presence  of 
unknown  AR  noise. 
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Figure  4.2.  Angle  estimates  obtained  from  the  experimental  data,  corresponding 
to  64  observation  intervals.  The  solid  lines  denote  the  means  and  the  dashed  lines 
denote  the  means  plus  and  minus  the  standard  deviations  of  the  angle  estimates.  The 
true  incident  angles  are  believed  to  be  61  =  —33°  and  62  =  —36°.  (a)  RELAX,  (b) 
CSM-ESPRIT. 


CHAPTER  5 
ANGLE  AND  POLARIZATION  ESTIMATION  WITH  A  COLD  ARRAY 


5.1     Introduction 

This  chapter  studies  the  advantages  of  an  arbitrary  linear  array  that  consists 
of  Co-centered  Orthogonal  Loop  and  Dipole  (COLD)  pairs.  By  using  the  COLD 
array,  the  performance  of  both  angle  and  polarization  estimation  can  be  significantly 
improved  as  compared  to  using  a  co-centered  crossed  dipole  (CCD)  array. 

The  case  where  all  incident  narrowband  electromagnetic  (EM)  plane  waves  are 
completely  polarized  is  considered.  A  completely  polarized  EM  wave  is  a  limiting  case 
of  a  more  general  type  of  EM  wave,  viz.  a  partially  polarized  EM  wave.  The  state  of 
polarization  of  a  partially  polarized  EM  wave  is  a  function  of  time  while  a  completely 
polarized  wave  has  a  fixed  state  of  polarization  (see  [102]  and  the  references  therein). 

We  present  an  asymptotically  statistically  efficient  signal  subspace-based  MODE 
algorithm  [25,  26]  for  both  angle  and  polarization  estimation.  Since  the  MODE  algo- 
rithm is  a  signal  subspace-based  approach,  it  is  asymptotically  statistically  efficient 
for  both  correlated  (including  coherent)  and  uncorrelated  incident  signals.  We  show 
with  numerical  examples  that  the  estimation  performance  of  MODE  is  better,  espe- 
cially for  highly  correlated  or  coherent  signals,  than  that  of  MUSIC  and  NSF  (noise 
subspace  fitting)  [65].  (We  remark  that  the  signal  subspace  eigenvector  based  MODE 
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algorithm  and  the  NSF  algorithm  are  asymptotically  statistically  equivalent  whenever 
the  signals  are  non-coherent.  For  coherent  signals,  MODE  remains  asymptotically 
statistically  efficient,  whereas  NSF  is  no  longer  asymptotically  statistically  efficient 
[76].  This  observation  suggests  that  when  the  correlation  coefficient  is  close/very 
close  to  one,  the  NSF  may  need  a  much  larger  number  of  data  samples  than  MODE 
to  converge  to  the  asymptotics,  and  hence  for  a  given  finite  N,  MODE  is  likely  to 
perform  better  than  NSF  in  such  a  case  of  highly  correlated  signals.) 

5.2     COLD  Array  and  Problem  Formulation 

Consider  a  2L-element  Hnear  array  consisting  of  L  Co-centered  Orthogonal  Loop 
and  Dipole  (COLD)  pairs  as  shown  in  Figure  5.1.  The  signal  received  from  each 
antenna  sensor  is  to  be  processed  separately  for  direction  and  polarization  estimation. 
The  /th  COLD  pair,  /  =  1, 2,  •  •  • ,  L,  has  its  center  on  the  j/-axis  at  an  arbitrary  y  =  6i. 
For  the  /th  COLD  pair,  the  dipole  parallel  to  the  z-axis  is  referred  to  as  the  z-axis 
dipole  and  the  loop  parallel  to  the  x-y  plane  as  the  x  -  y  plane  loop. 

Assume  K  (with  K  <  L)  narrowband  plane  waves  impinge  on  the  array  from 
angular  directions  described  by  9  and  </>,  where  9  and  4>  denote  the  azimuth  and 
elevation  angles,  respectively,  as  shown  in  Figure  5.1.  Furthermore,  suppose  each 
signal  is  a  completely  polarized  transverse  electromagnetic  wave  with  an  arbitrary 
eUiptical  electromagnetic  polarization  [103].  Assume  that  the  electric  field  of  an 
incoming  signal  has  transverse  components 


E  =  Eeee  +  ^^e^,  (5.1) 
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where  unit  vectors  e^,  e^,  and  -e^,  in  that  order,  form  a  right-hand  coordinate  system 
for  the  incoming  signals  and  £$  and  E^  are  the  horizontal  and  vertical  components  of 
the  electric  field  respectively.  In  general,  as  time  progresses.  Eg  and  E^  will  describe 
a  polarization  ellipse.  For  a  given  signal  polarization,  specified  by  constants  7  and  77, 
the  electric  field  components  are  given  by  (aside  from  a  common  narrowband  phase 
factor  So{t)  ^) 

Ee    =    £cos7,  (5.2) 

E^    =    Esin7e^^  (5.3) 

where  E  denotes  the  amplitude  of  the  incident  signal.  The  7  and  t/  can  be  used  to 
compute  a  and  /?,  which  are  the  ellipticity  and  orientation  angles  of  the  polarization 
ellipse,  respectively.  7  is  always  in  the  range  0  <  7  <  7r/2  and  rj  is  in  the  range 
— TT  <  7;  <  TT.  a  and  (3  can  also  be  used  to  compute  7  and  77  [31,  40]. 

Assume  that  each  dipole  in  the  array  is  a  short  dipole  (i.e.,  the  length  of  the 
dipole  is  equal  to  or  less  than  one-tenth  of  a  wavelength)  with  the  same  length  Lsd 
and  each  loop  is  a  small  loop  (i.e.,  the  perimeter  of  the  loop  is  equal  to  or  less  than 
three-tenths  of  a  wavelength)  with  the  same  area  Asi-  Thus  the  output  voltages  from 
each  dipole  and  loop  are  proportional  to  the  electric  field  components  parallel  to 
dipole  and  loop,  respectively. 


^For  a  narrowband  BPSK  (binary  phase-shift  keyed)  signal,  for  example,  «<,(<)  =  e''I'^°*+'''(*)l, 
where  Wo  is  the  carrier  frequency  and  V'(0  is  the  modulating  phase. 
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An  incoming  signal  described  by  arbitrary  electric  field  components  Eg  and  E^ 


can  be  written  as 


E  =  E 


(cos  7)69  +  {s'm^e^'')e^ 


(5.4) 


Let  us  define  the  spatial  phase  factor 


.2iv6,     .     , 
7— r — ^  smo 


(5.5) 


where  Aq  is  the  wavelength  of  the  signal.   The  effective  heights  of  the  short  dipoles 
and  small  loops  are  given  by  [104] 


/isd  =  isd  sin  ^, 


(5.6) 


and 


hsi  =  -J—; — smcp, 

^0 


(5.7) 


respectively.  Including  the  time  and  space  phase  factors  in  (5.4),  we  find  that  an 
incoming  signal  characterized  by  {9,  ^,7, 77,  E)  produces  a  signal  vector  in  the  COLD 
pair  centered  ai  y  =  61  as  follows: 


Mt) 


Xl{t) 
Xd{t) 


=  IJ-coLDEsoit)qi, 


(5.8) 
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where 


t^COLD 


j  ^^^  sin  (f)  cos  7 


-Lsd  s\n(f>smje 


jri 


(5.9) 


An  advantage  of  the  COLD  array  is  that  its  antenna  elements  are  not  sensitive  to 
the  azimuth  angle  6  of  the  signal  because  both  the  loops  and  dipoles  have  the  same 
sincj)  field  pattern,  as  may  be  seen  from  Equation  (5.9).  Hence  the  incoming  signal 
described  by  (5.4)  is  independent  of  6.  We  assume  that  the  antennas  and  the  incident 
signals  are  co-planar,  i.e.,  4>  =  90°.  Thus  (5.9)  becomes 


t^COLD 


Ve  cos  7 
Vrj,  sin  76-''' 


(5.10) 


where 


Ve 


.2irA, 


J- 


Ao    ' 
V^    =    —Lsd- 


(5.11) 
(5.12) 


Note  that  Ve  and  V^  represent  the  complex  voltages  induced  at  the  loop  and  dipole 
outputs  by  a  signal  with  a  unit  electric  field  parallel  to  the  loops  and  dipoles,  respec- 
tively. Let  s{t)  =  Eso{t)Ve  cos 'y.  The  Zi{t)  in  (5.8)  can  be  rewritten  as 


z{t)  =  us{t)qi, 


(5.13) 
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where 


u 


-^  tan  76^" 


(5.14) 


Assume  that  K  signals,  specified  by  incident  angles  ^i,  k  =  1,2,  •■•,A',  are 
incident  on  the  array.  In  addition,  we  assume  a  thermal  noise  voltage  vector  n/(f)  is 
present  at  each  output  vector  Z/(f).  The  n;(f)  are  assumed  to  be  zero-mean  circularly 
symmetric  complex  Gaussian  random  processes  that  are  statistically  independent  of 
each  other  and  to  have  covariance  matrix  cr^I,  where  I  denotes  the  identity  matrix. 

Under  these  assumptions,  the  total  output  vector  received  by  the  COLD  pair 
centered  at  y  =  ^/  is  given  by 


Z/(i)  =  ^Ufc5fc(i)g/fc  +  ni(0,      /  =  1,2,  •••,/>,  (5.15) 


where  u^  and  qik  are  given  by  (6.9)  and  (5.5),  respectively,  with  subscript  k  added  to 
each  angular  quantity.  Further,  Sk{t)  =  EkSok{t)Ve  cos 'yk,  where  EkSok{t)  denotes  the 
kth  narrowband  signal.  The  incident  signals  may  or  may  not  be  correlated  (including 
completely  correlated,  i.e.,  coherent)  with  each  other. 
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Let  z{t),  s{t),  and  n(^)  be  column  vectors  containing  the  received  signals,  inci- 
dent signals,  and  noise,  respectively,  i.e., 


z{t)  = 


z,{t) 
z,{t) 


,      s(0  = 


SK{t) 


,     n(i)  = 


Mi) 

Mt) 


The  received  signal  vector  has  the  form 


(5.16) 


z(^)  =  As{t)  +  nit), 


(5.17) 


where  A  is  a  2L  x  K  matrix 


A    = 


ai    R2    •  •  •    ^K 


=    [A®I]U, 


(5.18) 


with  (8)  representing  the  Kronecker  product, 


qu    qi2    •••    qiK 


921        922       •  •  •       q2K 


qu    qL2   •••   qLK 


(5.19) 
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and 


Ui 


0 


U 


(5.20) 


0 


ua 


Assume  that  the  element  signals  are  sampled  at  N  distinct  times  ^„,  n  =  1, 2,  •  •  • ,  A'^. 
The  random  noise  vectors  n(^„)  at  different  sample  times  are  assumed  to  be  in- 
dependent of  each  other.  The  problem  of  interest  herein  is  to  determine  the  az- 
imuth arrival  angles  Ok  and  the  states  of  polarization  described  by  (7^,  rik)  or  [ak.^k), 
fc  =  1, 2,  ■  •  • ,  /i',  from  the  measurements  z(f„),  n  =  1, 2,  •  •  • ,  A^. 

5.3     Angle  and  Polarization  Estimation  using  MODE 

The  MODE  [25,  26]  and,  in  a  related  form,  the  WSF  [23]  algorithms  were  derived 
for  angle  estimation  with  uniformly  polarized  arrays.  We  present  below  how  to  use 
the  signal  subspace-based  MODE  algorithm  with  the  COLD  array  for  both  angle  and 
polarization  estimation. 

Let 


N 


R=^E==(^n)z''(0, 


(5.21) 


n=\ 


where  (■)^  denotes  the  complex  conjugate  transpose  and  R  denotes  the  estimate  of 
the  following  array  covariance  matrix: 


Ti=E[z{t)z^{t)]. 


(5.22) 
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It  has  been  shown  in  [25,  23,  76]  that  an  asymptotically  (for  large  N)  statistically 

T 

efficient  estimator  of  the  angles  0  =  [  5)j ^    ^2,    •••,    ^a'  ]     and  the  polarization  pa- 


rameters r  = 


function: 


n,   ^2, 


rK 


can  be  obtained  by  minimizing  the  following 


/(0,r)  =  Tr 


.',  2 


-1, 


P^E.AA    Ef 


(5.23) 


-i-H 


where  the  symbol  P^  stands  for  the  orthogonal  projector  onto  the  null  space  of  A 
and  the  columns  in  Ej  are  the  signal  subspace  eigenvectors  of  R  that  correspond  to 
the  K  largest  eigenvalues  of  R,  with  K  =  min[A'^,  rank(S)].    Here  S  is  the  source 


covariance  matrix, 


S  =  E{s{tn)s"{tr,)}. 


(5.24) 


Assume  that  K  is  known.  (If  K  is  unknown,  it  can  be  estimated  from  the  data  as 
described,  for  example,  in  [43].)  Note  that  if  no  components  of  the  signal  vector  s{t) 
are  fully  correlated  to  one  another,  then  K  =  K  (provided  A'^  >  A').  Further,  the  A 
in  (5.23)  is  a  diagonal  matrix  with  diagonal  elements  Ai  >  A2  >  •  •  •  >  A^.,  which  are 
the  K  largest  eigenvalues  of  R,  and 


A  =  A  -  a% 


(5.25) 


where 


a'  = 


1 


2L 


2L  -  A' 


>  E  A» 


1 


=A-+1 


2A  -  A 


tr(R)-£X- 


(5.26) 
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We  show  below  that  we  can  concentrate  out  r  first  and  hence  reduce  the  dimen- 
sion of  the  parameter  space  over  which  we  need  to  search  to  minimize  (5.23). 
It  has  been  shown  (in  Appendix  B)  that 


>^  =  I-Pa«i  +  P    .. 


(Al«®I)V' 


(5.27) 


where 


At  =  (A^A)-'a^, 


(5.28) 


and 


V  = 


vi  0 


0  VA- 


(5.29) 


with 


Vfc  = 


,     k  =  l,2,---,K. 


Thus  minimizing  f{9,  r)  in  (5.23)  is  equivalent  to  minimizing 


(5.30) 


/(0,r)  =  Tr 


-P-«^  +  P.aKi,vP^^^"^^ 


-1^ 


(5.31) 


Let 


W=  |v^[(A^A)  '(8)I]v|     , 


(5.32) 
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be  formed  from  some  consistent  estimates  of  6  and  r.  Since  P     i     ^.^  E^  =  0{l/\/N), 

(A  I  "giljV 

|V^[(A^A)"^  ®  I]V|"    can  be  replaced  by  W  without  affecting  the  asymptotics  of 
the  MODE  estimator  [25,  26].  Then  we  have 


f{e,r)  =  h{e)  +  f2ie,r), 


(5.33) 


where 


fi{B)  =  -Tr 


(A®I){(A^A)  '®l}(A^®I)E,A  A  ^Ef 


(5.34) 


and 


/2(^,r)  =  Tr 


^  2 

(At^  ®  I)VWV^(At  ®  i)E^A  A~'Ef 


(5.35) 


The  MODE  estimates  {^,  r}  are  obtained  by  minimizing  /,  i.e., 


{^,r}  =  argrnin[/a(0)  +  /2(0,r)] 

tf,T 


(5.36) 


Let 


A/j  —  odd  columns  of  (A^^  ®  I), 


(5.37) 


and 


Ai,  =  even  columns  of  (A'^  ®  I). 


(5.38) 


Let  V/i  and  V„  be  the  following  K  x  K  diagonal  matrices: 
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V;,  =  diag|    ^rl,    ■ 


'  K 


(5.39) 


and 


V„  =  I. 


(5.40) 


Then 


(At^  ®  I)V  =  A,,V,,  +  A,V,. 


(5.41) 


Thus  /2(0,r)  in  (5.35)  can  be  rewritten  as 


,'■,  2 


+Tr  V^AfE,A'A"'Ef A,V,W 
+Tr  [vf  Af  E,A'A"'Ef  A,V,W 
+Tr  V,^AfE3A'A"'Ef A,V,W 


(5.42) 


Since  Vh  and  V„  are  diagonal  matrices,  Equation  (5.42)  can  be  written  in  the  fol- 
lowing matrix  form: 


f2{0,r) 


r 

-1 

V/^ 

vf 

e-^ 

Q(^) 

'- 

-^ 

e 

(5.43) 
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where 


Q{0) 


.',  2 


(A^E.A  A"'Ef  A,)  0  W^    (Af  E,A  A'^Ef  A,)  ©  W^ 
(Af  E,rA"'Ef  A,)  0  W^    (Af  E.A'A^'Ef  A,)  0  W^ 


(5.44) 


with  0  denoting  the  Hadmard-Schur  matrix  product  (i.e.,  the  elementwise  multipli- 
cation), 


Vh  = 


-\T 


-rT 


and 


iT 


(5.45) 


1     ■••     1 


(5.46) 


Note  that  the  polarization  parameters  are  contained  only  in  v/j.  By  setting  8/2/ dvh  = 


0,  we  obtain 


v,  =  -q^\e)Q2{0)e. 


(5.47) 


Using  (5.47)  in  (5.43)  gives 


f3{0)  =  e'  Q3{e)  -  Q^  {e)Q-\e)Q2{0) 


(5.48) 


which  is  a  concentrated  function  depending  only  on  0.  The  MODE  estimates  {^,f} 

can  be  obtained  by 

0  =  ^vgmm[fr{0)  +  fs{0)],  (5.49) 

0 
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and  using  9  in  (5.47)  to  obtain  r. 

To  summarize,  we  have  the  following  MODE  algorithm  for  angle  and  polariza- 
tion estimation: 

Step  1:  Obtain  initial  estimates  of  6  and  r  (see  the  discussions  below). 

Step  2:  Determine  8  by  minimizing  fi{d)  +  h{0)  as  shown  in  Equation  (5.49) 
with  W  in  (5.32)  formed  from  the  initial  estimates  obtained  in  Step  1. 

Step  3:  Calculate  r  by  using  the  6  obtained  in  Step  2  in  (5.47). 

Step  4:  Determine  the  7  and  r)  from  r  with 


7fc    =    tan  ^ 


TkVe 


r]k    =     -arg     -- 


(5.50) 


,     A;  =  1,2,---,A'.  (5.51) 


For  signals  that  are  not  highly  correlated  or  coherent  with  each  other,  the  initial 
estimates  of  0  and  r  in  Step  1  may  be  obtained  by  using  MUSIC  [84],  which  requires 
a  one-dimensional  search  over  the  parameter  space.  For  highly  correlated  or  coherent 
signals,  the  initial  estimate  of  6  may  be  determined  by  setting  W  =  I  and  minimizing 
/i(^)  +  /3(^)  as  shown  in  Equation  (5.49).  The  initial  estimate  of  r  can  be  calculated 
by  using  the  initial  estimate  of  9  in  (5.47).  The  initial  estimates  obtained  by  using 
MUSIC  for  non-coherent  signals  or  MODE  with  W  =  I  are  known  to  be  consistent 
[26,  71]. 
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5.4     Statistical  Performance  Analysis 

We  present  below  the  asymptotic  (for  large  N)  statistical  performance  of  MODE 
for  both  direction  and  polarization  estimation  with  the  COLD  array. 

Before  we  present  the  analysis  results,  however,  we  first  describe  the  method  we 
use  to  describe  the  accuracy  of  the  polarization  estimates.  For  reasons  discussed  in 
[31],  we  define  the  polarization  estimation  error  to  be  the  spherical  distance  between 
the  two  points  M  and  M  on  the  Poincare  sphere  that  represent  the  actual  state  of 
polarization  (7,7;)  and  the  estimated  state  of  polarization  (7,17),  respectively.  Let  ( 
be  the  angular  distance  between  M  and  M.  Then  [31] 

cos  C  =  cos  27  cos  27  +  sin 27  sin 27  cos(7/  -  rj).  (5.52) 

where  (  is  always  in  the  range  0  <  (  <  ^r. 

Applying  the  first-order  approximation  to  the  left  side  of  5.52  yields 

Cl  =  4(7.  -  Ik)'  +  ^m\2^k){flk  -  m)'-  (5.53) 

The  asymptotic  variances  of  the  polarization  estimates  are  obtained  with  (5.53)  and 
the  accuracy  results  on  7  and  t)  given  below. 
Let 


0^    7^    r7^ 


T 


(5.54) 
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It  follows  from  [25,  76]  that  the  asymptotic  (for  large  A^)  statistical  distribution  of  f 
is  Gaussian  with  mean  r  and  covariance  matrix  equal  to  the  corresponding  stochastic 
Cramer- Rao  bound  (CRB),  CRB.  The  ijth  element  of  CRB"^  is  given  by 


CRB-^]^.^.  =  ^Re  [tr{Af  P^A.-SA^'r-iAS}]  ,  (5.55) 


where  Ai  =  dA/dri  with  r,  denoting  the  ith.  element  of  r. 

5.5     Numerical  Results 

We  present  below  several  examples  showing  the  performance  of  using  the  MODE 
algorithm  with  the  COLD  array  and  comparing  the  asymptotic  statistical  perfor- 
mance analysis  results  with  the  Monte-Carlo  simulation  results.  We  compare  MODE 
with  MUSIC  and  NSF  for  both  angle  and  polarization  estimation.  The  simulation 
results  were  obtained  by  using  50  Monte-Carlo  simulations.  In  the  examples,  we  as- 
sume that  there  are  K  =  2  incident  signals  and  both  signals  are  assumed  to  have  the 
same  amplitude  Ek,  such  that  \VeEk\  =  \V^Ek\  =  I,  k  =  1,2.  Hence,  the  signal-to- 
noise  ratio  (SNR)  used  in  the  simulations  is  -lOlogjocr^  dB.  The  array  is  assumed 
to  have  L  =  8  COLD  pairs  that  are  uniformly  spaced  with  the  spacing  between  two 
adjacent  COLD  pairs  equal  to  a  half  wavelength.  We  also  compare  the  estimation 
performance  of  using  the  COLD  array  with  that  of  using  a  CCD  array  with  the 
same  array  geometry.  The  CCD  array  consists  of  crossed  y-  and  z-a,xes  dipoles.  The 
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counterpart  of  Equation  (5.9)  for  the  CCD  array  can  be  written  as 


t^CCD 


Lsd  cos  7  cos  d  +  Lsd  sin  76-'''  sin  6  cos  </> 
—Lsd  sin  (?!)  sin  76'''' 


(5.56) 


In  the  following  examples,  the  antennas  and  the  incident  signals  are  assumed  to  be 
co-planar,  i.e.,  (b  =  90°,  for  both  the  COLD  array  and  CCD  array.    For  (j)  =  90°, 


(5.56)  becomes 


l^CCD  — 


Lsd  cos  7  cos  6 
-Lsdsm'ye^'' 


(5.57) 


(We  remark  that  if  the  antennas  and  the  incident  signals  are  not  co-planar,  we  will 
need  two-dimensional  CCD  or  COLD  arrays  for  angle  and  polarization  estimation, 
which  is  the  case  not  considered  herein.  For  this  case,  however,  the  COLD  array  will 
not  always  perform  better  than  the  CCD  array.) 

First,  we  present  two  examples  that  illustrate  how  the  angle  separation  between 
the  two  incident  signals  affects  both  the  direction  and  polarization  estimates.  We 
begin  with  the  case  of  two  signals  with  identical  circular  polarizations  (qi  =  02  = 
45°).  Figure  5.2  shows  the  root-mean-squared  errors  (RMSEs)  of  the  estimates  of 
the  first  signal  as  a  function  of  angle  separation  A^  when  two  correlated  signals 
with  correlation  coefficient  0.99  arrive  at  the  array  from  angles  9i  =  -A6/2  and 
62  =  Ae/2.  We  note  that  MODE  performs  better  than  MUSIC  and  NSF.  Further, 
MODE  achieves  the  best  possible  unbiased  performance,  i.e.,  the  corresponding  CRB, 
as  the  angle  separation  increases.    Because  the  signals  arrive  from  angles  near  the 
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broadside  of  the  arrays,  the  CRBs  for  the  COLD  and  CCD  arrays  are  similar.  This 
case  corresponds  to  small  incident  angles,  which  make  HccD  similar  to  fJ-coLD  ^ 
may  be  seen  by  comparing  Equations  (6.5)  and  (5.57). 

In  Figure  5.3,  we  consider  the  case  where  the  signals  with  identical  horizontal 
polarizations  {a,  =02  =  0°,  l3i  =  ^2  =  0°)  arrive  from  angles  away  from  the  broadside 
of  the  array.  In  this  case,  the  CCD  array  is  outperformed  by  the  COLD  array. 
This  result  occurs  because  the  signal  outputs  at  the  j/-axis  dipoles  are  attenuated 
by  a  factor  of  cos  ^  for  the  CCD  array  (see  Equation  (5.56)).  For  the  COLD  array, 
however,  the  signal  outputs  at  both  the  dipoles  and  the  loops  are  independent  of  the 
incident  angle  9  (see  (6.5)).  Note  also  that  the  RMSEs  of  the  angle  estimates  first 
decrease  and  then  increase  even  as  the  angle  separation  increases.  This  result  occurs 
because  the  incident  angle  of  the  second  signal  approaches  90°  for  very  large  A^  and 
the  RMSEs  of  angle  estimates  are  approximately  proportional  to  —^  for  large  A^ 
[71].  We  note  again  that  MODE  gives  better  performance  than  MUSIC  and  NSF  and 
achieves  the  CRB  as  A^  increases.  We  have  also  found  that  the  CRBs  for  COLD 
array  are  much  lower  than  those  for  CCD  array,  especially  when  9  approaches  90°. 

Second,  we  consider  how  the  polarization  separation  affects  the  estimator  per- 
formance. Consider  the  case  where  two  incident  signals  with  a  correlation  coefficient 
0.99  arrive  from  angles  9i  =  50°  and  92  =  70°.  We  assume  that  the  corresponding 
ellipticity  angles  are  Qi  =  45°  -  Aa  and  02  =  45°  and  the  orientation  angles  are 
/3j  =  ^2  =  0°.  The  polarization  separation  between  the  two  polarization  states  is 
2Aq.  Figure  5.4  shows  the  RMSEs  of  the  direction  and  polarization  estimates  as  a 
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function  of  Aa  for  the  second  signal.  Note  that  MODE  performs  better  than  MUSIC 
and  NSF  since  the  signals  are  highly  correlated.  Since  the  second  signal  arrives  from 
a  large  angle  away  from  the  broadside  of  the  array,  the  performance  of  the  COLD 
array  is  again  much  better  than  that  of  the  CCD  array. 

Third,  we  consider  how  the  correlation  coefficient  between  the  two  incident 
signals  affects  both  the  direction  and  polarization  estimates.  Figure  5.5  shows  the 
RMSEs  of  the  direction  and  polarization  estimates  as  a  function  of  the  correlation 
coefficient.  The  two  signals  arrive  from  611  =  -6°  and  02  =  6°  and  we  have  ai  =  a2  = 
45°  and  0i  =  1^2  =  0°.  We  note  that  for  highly  correlated  or  coherent  signals,  the 
performance  of  MODE  is  much  better  than  that  of  MUSIC  and  NSF. 

Fourth,  we  consider  the  performance  of  the  estimators  as  a  function  of  SNR. 
In  Figure  5.6,  we  consider  the  case  where  6^  =  50°,  02  =  70°,  Oi  =  0:2  =  0°, 
^1  =  /32  =  0°,  correlation  coefficient  =  0.99,  and  A^  =  400.  In  this  case,  the  advantage 
of  MODE  over  MUSIC  and  NSF  is  increasingly  evident  at  low  SNR.  We  note  again 
that  the  COLD  array  gives  much  lower  CRB  than  the  CCD  array,  especially  for  angle 
estimation. 

Fifth,  we  consider  the  performance  of  the  MODE  estimator  as  a  function  of  the 
number  of  snapshots  N.  We  consider  the  case  where  9i  =  50°,  $2  -  70°,  (y.x=  0L2  — 
45°,  ^1  =  /32  =  0°,  correlation  coefficient  =  0.99.  and  SNR  =  10  dB.  It  can  be  seen 
from  Figure  5.7  that  although  MODE  is  an  asymptotically  (for  large  N)  statistically 
efficient  estimator,  MODE  estimates  achieve  the  CRB  even  for  very  moderate  N. 
MODE  again  performs  better  than  MUSIC  and  NSF. 
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Finally,  we  remark  that  although  we  make  the  Gaussian  noise  assumption  in  the 
problem  formulation,  the  Gaussian  noise  assumption  is  not  critical  to  the  algorithms. 
As  a  matter  of  fact,  the  asymptotic  distributional  properties  of  all  of  the  algorithms 
are  most  likely  to  remain  the  same  even  in  the  non-Gaussian  noise  case.  However, 
of  course,  in  the  latter  case,  MODE  is  no  longer  asymptotically  statistically  efficient 
and  the  CRB  given  in  (5.55)  is  no  longer  the  true  CRB  matrix,  but  rather  the 
asymptotic  theoretical  performance  (TP)  of  MODE.  Figure  5.8  shows  an  example 
when  the  Gaussian  noise  assumption  is  violated.  Figure  5.8  is  the  same  as  Figure 
5.7  except  that  the  noise  is  the  contaminated  Gaussian  noise  [105].  More  specifically, 
the  probability  density  function  of  n{t)  has  the  form  (1  -  e);V(0,  a^T)  +  ejV(0,  9aH), 
where  A/'(a,B)  denotes  the  Gaussian  probability  density  function  with  mean  a  and 
covariance  matrix  B.  The  SNR  here  is  defined  as  -101ogiok^(l  +  8e)]  dB.  We 
used  e  =  0.7  in  the  example.  Note  that  MODE  can  also  achieve  its  asymptotic 
theoretical  performance  for  moderately  large  N  for  this  case  and  MODE  has  the  best 
performance. 

5.6     Conclusions 

We  have  presented  a  Co-centered  Orthogonal  Loop  and  Dipole  (COLD)  array. 
We  have  shown  that  with  the  COLD  array,  both  angle  and  polarization  estimation 
performance  can  be  significantly  improved  as  compared  to  using  a  similar  crossed 
dipole  array.  We  have  also  described  how  to  use  the  asymptotically  statistically  ef- 
ficient MODE  algorithm  for  both  angle  and  polarization  estimation  with  the  COLD 
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array.  Numerical  examples  have  been  given  to  show  that  MODE  gives  more  accu- 
rate angle  and  polarization  estimates  than  MUSIC  and  NSF,  especially  for  highly 
correlated  or  coherent  signals. 
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Figure  5.1.  A  linear  COLD  array. 
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Figure  5.2.  Root-mean-squared  errors  (RMSEs)  of  estimates  versus  A^  for  the  first 
of  the  two  signals  when  ^i  =  -A9/2,  02  =  -FA^/2,  a^  =  a2  =  45°,  /3i  =  ^2  =  0°, 
correlation  coefficient  =  0.99,  N  =  400,  and  SNR  =  10  dB.  (The  CRBs  for  the  CCD 
array  nearly  coincide  with  those  for  the  COLD  array.)  (a)  Direction  estimates,  (b) 
Polarization  estimates. 
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Figure  5.3.  Root-mean-squared  errors  (RMSEs)  of  estimates  versus  A^  for  the  second 
of  the  two  signals  when  Qx  =  50°,  02  =  50°  +  A^,  a^  =  02  =  0\  ^i  =  ^2  =  0°, 
correlation  coefficient  =  0.99,  A'^  =  400,  and  SNR  =  10  dB.  (a)  Direction  estimates, 
(b)  Polarization  estimates. 
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Figure  5.4.  Root-mean-squared  errors  (RMSEs)  of  estimates  versus  Aa  for  the  second 
of  the  two  signals  when  9i  ==  50° ,  ^2  =  70° ,  ai  =  45°  -  Aa  and  aa  =  45° ,  ^^i  =  /?2  =  0° , 
correlation  coefficient  =  0.99,  N  =  400,  and  SNR  —  10  dB.  (a)  Direction  estimates, 
(b)  Polarization  estimates. 
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Figure  5.5.  Root-mean-squared  errors  (RMSEs)  of  estimates  versus  source  correlation 
coefficient  for  the  first  of  the  two  signals  when  9i  =  —6°,  ^2  —  6°,  oi  =  a2  =  45°, 
^1  =  ^2  =  0°,  N  =  400,  and  SNR  =  10  dB.  (The  CRBs  for  the  CCD  array  nearly 
coincide  with  those  for  the  COLD  array.)  (a)  Direction  estimates,  (b)  Polarization 
estimates. 
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Figure  5.6.  Root-mean-squared  errors  (RMSEs)  of  estimates  versus  SNR  for  the 
second  of  the  two  signals  when  Ox  =  50°,  82  =  70°,  ai  =  a2  =  0°,  ^1  =  ^2  =  0°, 
correlation  coefficient  =  0.99,  and  N  =  400.  (a)  Direction  estimates,  (b)  Polarization 
estimates. 
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Figure  5.7.  Root-mean-squared  errors  (RMSEs)  of  estimates  versus  A''  for  the  second 
of  the  two  signals  when  Oi  ==  50°  and  6*2  =  70°,  qi  =  ^2  =  45°,  /9i  =  /?2  =  0°,  correla- 
tion coefficient  =  0.99,  and  SNR  =  10  dB.  (a)  Direction  estimates,  (b)  Polarization 
estimates. 
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Figure  5.8.  Root-mean-squared  errors  (RMSEs)  of  estimates  versus  N  for  the  second 
of  the  two  signals  in  the  presence  of  contaminated  Gaussian  noise  when  6i  =  50°  and 
92  =  70°,  Qi  =  Q2  =  45°,  I3i=  132  =  0°,  correlation  coefficient  =  0.99,  and  SNR  =  10 
dB.  (a)  Direction  estimates,  (b)  Polarization  estimates. 


CHAPTER  6 
PARAMETER  ESTIMATION  USING  RELAX  WITH  A  COLD  ARRAY 


6.1     Introduction 

The  RELAX  algorithm  described  in  [45]  is  both  conceptually  and  computation- 
ally simple:  its  implementation  mainly  requires  a  sequence  of  fast  Fourier  transforms 
(FFTs).  However,  the  RELAX  algorithm  has  been  used  only  for  either  temporal  spec- 
tral analysis  [45]  or  angle  and  waveform  estimation  with  a  uniform  Hnear  array  [35]. 
In  this  Chapter,  we  derive  a  RELAX  approach  for  angle,  polarization  and  waveform 
estimation  of  completely  polarized  narrowband  electromagnetic  plane  waves  arriving 
at  a  polarization  sensitive  array,  namely  a  uniform  linear  COLD  (co-centered  orthog- 
onal loop  and  dipole)  array  as  described  in  Chapter  5.  This  new  RELAX  estimator 
also  requires  only  a  sequence  of  FFTs.  The  performance  of  this  RELAX  estima- 
tor will  be  compared  with  that  of  the  MODE  estimator  in  Chapter  5  via  numerical 
examples. 

6.2     Problem  Formulation 

The  COLD  array  we  consider  herein  is  a  uniform  linear  COLD  array  with  spacing 
8  between  the  adjacent  sensor  pairs,  as  shown  in  Figure  6.1.  For  a  given  signal 
impinging  on  the  array  from  angles  6  and  (/>,  where  9  and  (j)  denote  the  azimuth  and 
elevation  angles,  respectively,  with  its  state  of  polarization  specified  by  constants  7 
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and  T],  the  electric  field  components  are  given  by 


Eg      —      E  COS  7, 


Ea,    =    £  sin  76-''', 


(6.1) 
(6.2) 


where  E  denotes  the  amplitude  of  the  incident  signal.  The  7  and  rj  can  be  used 
to  compute  a  and  ^  [31,  40],  which  are  the  ellipticity  and  orientation  angles  of  the 
polarization  ellipse,  respectively.  An  incoming  signal  characterized  by  {9,(f),j,rj,E) 
produces  a  signal  vector  in  the  COLD  pair  centered  at  y  =  (/  -  1)^  as  follows: 


zi{t)  = 


x,{t) 

Xd{t) 


■  2w(l-l)6     .     „ 


where  Aq  is  the  wavelength  of  the  signal,  and 


(6.3) 


f^COLD  — 


j^^^sin(^cos7 
—  Lgd  sin  (fysiwye^^ 


(6.4) 


with  Lsd  and  Asi  denoting  the  length  of  the  dipole  and  the  area  of  the  loop,  respec- 


tively. 


For  the  case  that  the  antennas  and  the  incident  signals  are  co-planar,  i.e.,  (j) 


90°,  ((6.4)  becomes 


t^COLD 


Ve  cos  7 
V^  sin  76-''' 


(6.5) 
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where 


Vs 


Ao 
=    —Lsd- 


(6.6) 
(6.7) 


Note  that  Vg  and  V^  represent  the  complex  voltages  induced  at  the  loop  and  dipole 
outputs  by  a  signal  with  a  unit  electric  field  parallel  to  the  loops  and  dipoles,  respec- 
tively. Let  s{t)  =  Eso(t)Ve  cos -f.  The  zi{t)  in  (6.3)  can  be  rewritten  as 


.27r(l-l)«     .     a 


(6.8) 


where 


u 


^  tan  76^'^ 


Ve 


(6.9) 


Assume  that  K  signals,  specified  by  incident  angles  6k,  and  polarization  states  r^, 
k  —  1, 2,  •  •  • ,  A',  are  incident  on  the  COLD  array.  In  addition,  we  assume  a  thermal 
noise  voltage  vector  n;(^)  is  present  at  each  output  vector  z;(f).  The  {n;(i)};^^  are 
assumed  to  be  zero-mean  random  vectors  that  are  statistically  dependent  of  each 


other  and  to  have  an  unknown  covariance  matrix. 
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Under  these  assumptions,  the  total  output  vector  received  by  a  uniform  linear 
array  consisting  of  L  pairs  of  COLD  can  be  written  as: 


+  n(0,  (6.10) 


where  Zsi,k{t),  the  output  vector  from  the  L  small  loops  by  the  fcth  signal,  is  given  by 


z,;,fc(0  =  a{Ok)sk{t), 


(6.11) 


and  Zsd,k{t),  the  output  vector  from  the  L  small  dipoles  by  the  A;th  signal,  is  given  by 


Zsd,fc(0  =  rka.{6k)sk{t). 


(6.12) 


with  a(^fc),  k  =  1, 2,  •  ■  ■ ,  iC,  being  an  X  x  1  vector  for  the  kth  signal.  a(^A;)  has  the 
form 


a(^fc)  = 


e   -^0        * 


^,^(L-l)sine, 


fc  =  l,2,---,K  (6.13) 


n(t)  is  a  zero-mean  random  vector  with  an  unknown  covariance  matrix  Q.  We  define 
the  angle  and  polarization  vectors  as: 


-iT 


e  = 


r    = 


Or    02    ■■■    6k 


ri     r2     ■■■     rK 


(6.14) 
(6.15) 


Then  the  received  signal  vector  can  be  written  as: 
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z{t)  =  A(0,r)s{t)  +  n{t), 


(6.16) 


where  A{S,r)  is  the  2L  x  K  matrix 


A{e,v) 


i(^i,ri)    a{e2,r2)    ■■■    a{9K,rK) 


a(^i)        a(^2) 
ria(^i)    r2a(^2) 


a(eA-) 


and 


s{t)  = 


Slit) 

S2{t) 
SK{t) 


nit) 


Mt) 

n2{t) 

Mt) 


(6.17) 


(6.18) 


Assume  that  the  element  signals  are  sampled  at  N  distinct  times  t  =  1,2,-  ■  ■  ,N.  The 
random  noise  vectors  n{t)  at  different  sample  times  are  assumed  to  be  independent  of 
each  other.  The  problem  of  interest  herein  is  to  determine  the  azimuth  arrival  angles 
9k,  the  states  of  polarization  r^  (described  by  {^k,Vk)  or  (ajt,/9fc)),  and  the  signal 
waveforms  Sk{t),  k  =  1,2,-  ■  ■  ,K,  t  =  1,2,- ■  ■  ,N,  from  the  measurements  {z(t)}^i. 

6.3     Parameter  Estimation  with  RELAX 

The  RELAX  algorithm  has  been  shown  to  be  an  efficient  algorithm  for  temporal 
spectral  analysis  and  angle  and  waveform  estimation  with  a  uniform  linear  array.  We 
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describe  below  how  the  RELAX  algorithm  can  be  derived  for  parameter  estimation 
with  a  uniform  linear  COLD  array. 

The  RELAX  algorithm  minimizes  the  nonlinear  least-squares  (NLS)  criterion 
for  the  data  model  in  (6.16),  which  is  given  by 


N 

I 


Ci  [0,  r,  s(l),  s(2),  •  •  • ,  s{N)]  =  YlKt)  -  A(e,  r)s{t)f[^{t)  -  A(0,  r)s(i)].     (6.19) 


If  the  noise  n{t)  were  spatially  white  and  Gaussian,  minimizing  (6.19)  would  have 
yielded  the  deterministic  maximum  likelihood  estimates  of  the  incident  angles,  po- 
larization states,  and  signal  waveforms  [20]. 

The  NLS  criterion  for  estimating  the  kih  signal  parameters  can  be  written  as: 


A' 


C2[^fc,  r-fc,  Sfc(l),  5fc(2),  •  •  • ,  Sk{N)]  =  Yl  [zfc(0  -  a(^fc>  n)sk{t)f  [zk{t)  -  a(^,,  rk)sk{t)] , 

(6.20) 

where 

Zk{t)  =  zit)-    Yl    Bi{e,,r,)si{t\  (6.21) 

and       where       K       is       an       intermediate       number       of       signals,         and 
{Oi,  ri,  5i(l),  Si{2),  •  •  • ,  Si{N)}^^i^i^k  are  assumed  to  have  been  previously  estimated. 
Minimizing  C2  with  respect  tosfc(l),---,5fc(A^)  gives  their  estimates 


,  a"{ek,rk)zk{t) 

Sk{t)      = 


L{l  +  \rk\') 


ei"iek)zsi,k{t)  +  r^a^(gfc)z,,,,(0 


L{l  +  \rk\') 


,     t  =  l,2,..-,A^.(6.22) 

h=Sk,^k=fk 


The  estimates  of  dk  and  r*;,  A;  =  1, 2,  •  •  • ,  K,  are  obtained  by  minimizing 
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N 


C3{0k,rk)  =  Yl 


t=i 


I- 


i{ek,rk)^"{0k,rk) 


L{l  +  In?) 


Zk{t) 


(6.23) 


or  maximizmg 


1  +  tm'  ^  ' 


(6.24) 


which  can  be  written  as: 


CsiOk,  n)  =  :r^-i^[e/(^0  +  \rk\'ed{ek)  +  neu{e,)  +  r*£;,(^,)],  (6-25) 

1  +  VkV 


where 


N 


CiiOk) 

=    Yl  a^(0fc)z./,.(^)    , 
t=i 

(6.26) 

e40k) 

t=i 

(6.27) 

eid{Ok) 

=  )^{[^"{Ok)zsi,k{m^''i0k)zsdAi)r}- 

(6.28) 

t=l 


Setting  to  zero  the  derivatives  of  C5  with  respect  to  the  real  and  imaginary  parts  of 
Tfe  yields 

rleidiek)  +  rk[e,{ek)  -  €^{9^)]  -  e^^ie^)  =  0,  (6.29) 

which  gives  the  estimate  of  r^  of  r^  as 


rk  = 


ediOk)  -  ei{0k)  +  ^/jeijek)  -  e^jOkW  +  Meid{9k) 


'^£ld{9k) 


(6.30) 


6k=e 
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Note  that  we  chose  the  first  root  (plus  sign)  of  Equation  (6.29)  as  the  ft  shown  in 
(6.30)  because  only  this  root  guarantees  that  the  second  derivatives  of  C5  with  respect 
to  the  real  and  imaginary  parts  of  r^  are  negative.  In  other  words,  this  root  assures 
the  maximization  of  C5  with  respect  to  rk- 

Inserting  (6.30)  into  (6.25)  leads  to  a  function  CeiOk),  which  depends  on  ei{9k), 
SdiOk)  and  eid{Bk)  that  depend  only  on  6k-  Thus  the  estimates  Ok  of  Ok  can  be  obtained 

by 

^,  =  argmax{C6(^fc)}.  (6.31) 

Note  that  since  C^iOk)  is  a  combination  of  ei{Ok),  ed{Ok)  and  eid{Ok)-,  it  can  be 
efficiently  computed  by  using  FFT  with  the  data  sequences  Zsi,k{t)  and  Xsd,k{t)  padded 
with  zeros.  Hence  Ok  can  be  obtained  as  the  location  of  the  dominant  peak  of  Ce{Ok)- 
Once  Ok  is  determined,  7"'^  can  be  readily  calculated  from  (6.30).  Once  {^fe,^.}  is 
determined,  Sk{t)  can  be  readily  computed  by  using  (6.22). 

Assume  that  the  number  of  incident  signals  K  is  already  known  or  estimated. 
Then  the  RELAX  algorithm  has  the  following  steps: 

Step  (1):  Let  K  =  1.  Obtain  ^1,  ri  and  Si{t)  from  Zi(i)  as  in  (6.31),  (6.30)  and 
(6.22),  respectively,  t  =  l,2,---,N. 

Step  (2):  Let  K  =  2.  Compute  Z2(t)  with  (6.21)  by  using  Oi,  fk  and  Si{t)  obtained 
in  Step  (1),  t  =  l,2,---,iV.  Obtain  ^2,  h  and  S2{t)  from  Z2{t)  as  in  (6.31), 
(6.30)  and  (6.22),  respectively,  t  =  1,2,  ■  ■  ■  ,7V. 
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Next,  compute  Zi{t)  with  (6.21)  by  using  §2,  h  and  S2{t)  and  redetermine  6'i,  h 
and  51  (n)  from  Zi(i)  as  in  (6.31),  (6.30)  and  (6.22),  respectively,  i  =  1, 2,  •  •  • ,  iV. 
Iterate  the  previous  two  substeps  until  "practical  convergence"  is  achieved  (to 
be  discussed  later  on). 

Step     (3):        Let     K      ^      3.  Compute    z^(t)     with     (6.21)     by    using 

{6i,ri,Si[\),Si{2),-  ■  ■  ,Si{N)}}^^  obtained  in  Step  (2).  Obtain  ^3,  ^3  and  h{t) 
from  Z3(t)  as  in  (6.31),  (6.30)  and  (6.22),  respectively,  ^  =  1,2,  •  •  •  ,  A^. 

Next,  compute  Zi(i)  with  (6.21)  by  using  {^i,fi,5,(l),ii(2),  •  •  •  ,ii(A^)}f=2  and 
redetermine  ^1,  h  and  h{t)  from  Zi(t),  t  =  1,2,---,A^.   Then  compute  Z2(0 
with  (6.21)  by  using  {^i,f„Si(l),Si(2),  •  •  •  ,5i(iV)}i=i,3  and  redetermine  ^2,  r2 
and  S2{t)  from  Z2(i),  t  =  1,2,-  ■  ■  ,N. 
Iterate  the  previous  three  substeps  until  "practical  convergence"  is  achieved. 

Steps  (4)  -  (K):  Continue  similarly  until  K  is  equal  to  the  known  or  estimated 
number  of  signal  sources.  (See  [45]  for  a  possible  approach  to  estimate  K  when 
it  is  unknown.) 

Final  Step:  Determine  7  and  r)  from  r  with 

rkVe 


^k  =  tan  ^ 


14 


(6.32) 


and 


^fc^arg^^V      k  =  1,2,  ■■■,!<.  (6.33) 
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The  "practical  convergence"  in  the  iterated  substeps  of  the  above  RELAX  al- 
gorithm may  be  determined  by  checking  the  relative  change  of  the  cost  function  in 
(6.19)  between  two  consecutive  iterations.  In  our  numerical  examples,  to  be  presented 
in  the  next  section,  we  terminate  the  iterations  when  the  relative  change  is  less  than 
or  equal  to  e  =  10~^. 

6.4     Numerical  Results 

We  present  below  several  examples  showing  the  performance  of  the  RELAX 
algorithm  for  parameter  estimation  with  the  COLD  array.  We  also  compare  the 
performance  of  the  RELAX  estimator  with  that  of  the  MODE  estimator  described  in 
[34].  The  simulation  results  were  obtained  by  using  50  Monte-Carlo  simulations.  In 
the  examples,  we  assume  that  there  are  K  =  2  incident  signals  and  both  signals  are 
assumed  to  have  the  same  amplitude  Ek-,  such  that  IV^JStl  =  iV^^^fcl  =  1,  A:  =  1,2. 
The  noise  n;(i)  is  assumed  to  satisfy 

n;(0  =  -«in/-i(0  +  wKO,  (6-34) 

where  w;(i),  /  =  1, 2,  •  •  • ,  L,  are  independently  and  identically  distributed  zero-mean 
complex  Gaussian  random  vectors  with  covariance  matrix  cr^I,  with  I  denoting  the 
identity  matrix.  The  amplitudes  of  the  signal  waveforms  are  assumed  to  be  1.  The 
SNR  referred  to  below  is  defined  by 


SNR=101ogio(^-J^)    [dB].  (6.35) 
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The  array  is  assumed  to  have  L  =  8  COLD  pairs  that  are  uniformly  spaced 
with  the  spacing  between  two  adjacent  COLD  pairs  equal  to  a  half  wavelength.  The 
number  of  snapshots  taken  at  the  array  output  is  A^  =  2.  All  array  output  vectors 
are  zero  padded  to  2^^  before  using  FFT  in  the  RELAX  algorithm. 

The  polarization  estimation  error  is  defined  as  the  spherical  distance  between 
the  two  points  M  and  M  on  the  Poincare  sphere  that  represent  the  actual  state  of 
polarization  (7,7/)  and  the  estimated  state  of  polarization  (7,77),  respectively.  Let  ( 
be  the  angular  distance  between  M  and  M.  Then 

cos  (  =  cos  27  cos  27  +  sin  27  sin  27  cos(77  —  i}).  (6.36) 

where  (  is  always  in  the  range  0  <  C  <  vr. 

The  waveform  estimates  for  MODE  are  obtained  by  using 


s{t)  =  [A"{e,  f )  A(^,  r)]-^  a''(^,  r)z{t).  (6.37) 


The  root-mean-squared  error  (RMSE)  of  the  kth  waveform  estimate  is  defined  as 


The  first  example  illustrates  the  performance  of  the  estimators  as  a  function  of 
SNR  when  the  additive  noise  is  spatially  white  ,  i.e.,  ai  =  0  in  (6.34).  Consider  the 
case  where  the  two  signals  have  diffrent  linear  polarizations  (ofi  =  02  =  0°,/?i  =  0°, 
and  ^2  =  10°).  The  two  signals  are  assumed  to  be  correlated  with  correlation  coef- 
ficient 0.99  and  arrive  at  the  array  from  61  =  10°  and  62  =  22°.  Figures  6.2(a),  (b). 
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(c)  show  the  root-mean-squared  errors  (RMSEs)  of  angle,  polarization,  and  wave- 
form estimates,  respectively,  of  the  first  signal.  (The  results  for  the  second  signal 
are  similar.)  It  can  be  seen  from  Figure  1  that  RELAX  and  MODE  yield  similar 
performances  for  this  case  of  spatially  white  noise. 

The  second  example  considers  the  case  where  the  noise  is  spatially  colored  with 
Q^  =  -0.85e-'''/^  in  (6.34).  The  two  incident  signals  are  the  same  as  those  in  the  first 
example.  Figure  6.3  shows  the  RMSEs  of  the  parameter  estimates  of  the  first  signal 
as  a  function  of  SNR.  We  note  that  RELAX  outperforms  MODE  when  the  SNR  is 
low.  This  result  occurs  because  MODE  is  derived  under  the  assumption  of  spatially 
white  noise  while  the  true  noise  is  spatially  colored  and  unknown. 

Finally,  we  compare  the  computational  complexities  of  RELAX  and  MODE. 
The  amount  of  computations  needed  by  RELAX  depends  on  A'^,  K,  and  the  number 
of  zero-paddings.  RELAX  can  be  faster  than  MODE  for  small  N  and  K  and  at  a 
moderate  SNR  when  the  number  of  zero-paddings  required  is  small.  For  example, 
RELAX  requires  only  about  25%  of  the  amount  of  computations  required  by  MODE 
when  A^  =  1,  A'  =  2,  and  the  number  of  zero-paddings  is  2^^  When  A''  =  2,  RELAX 
requires  about  50%  of  the  amount  of  computations  needed  by  MODE.  However, 
MODE  can  require  much  fewer  computations  than  RELAX  for  large  A^  and  K  and 
at  a  high  SNR  when  the  required  number  of  zero-padding  is  large.  Even  so,  RELAX 
can  still  be  more  attractive  than  MODE  since  the  former  can  be  easily  implemented 
in  hardware  with  simple  FFT  chips  and  in  parallel. 
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6.5     Conclusions 

A  RELAX  algorithm  has  been  presented  for  angle,  polarization,  and  waveform 
estimation  of  narrowband  plane  waves  arriving  at  a  uniform  linear  COLD  array.  The 
RELAX  algorithm  only  requires  a  sequence  of  FFTs  and  is  robust  against  the  presence 
of  unknown  spatially  colored  noise.  We  have  also  presented  numerical  examples  to 
demonstrate  the  performance  of  RELAX  and  compared  it  with  that  of  MODE. 


Ill 


Figure  6.1.  A  uniform  linear  COLD  array. 
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X      RELAX 
O      MODE 


(a) 


(b) 


(c) 

Figure  6.2.  Root-mean-squared  errors  (RMSEs)  of  estimates  versus  SNR  for  the 
first  of  the  two  signals  in  the  presence  of  white  noise  when  6i  =  10°,  O2  =  22°, 
ai  =  02  =  0°,  (3i  =  0°,  ^2  =  10°,  correlation  coefficient  =  0.99,  and  N  =  2.  (a) 
Direction  estimates,  (b)  Polarization  estimates,  (c)  Waveform  estimates. 
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Figure  6.3.  Root-mean-squared  errors  (RMSEs)  of  estimates  versus  SNR  for  the  first 
of  the  two  signals  in  the  presence  of  unknown  AR  noise  when  0i  =  10°,  02  =  22°, 
ai  =  02  =  0°,  I3i  ^  0°,  /32  =  10°,  correlation  coefficient  =  0.99,  and  A^  =  2.  (a) 
Direction  estimates,  (b)  Polarization  estimates,  (c)  Waveform  estimates. 


CHAPTER  7 
CONCLUSIONS  AND  FUTURE  WORK 


7.1     Summary 

Most  existing  high  resolution  array  processing  algorithms  are  very  sensitive  to 
the  presence  of  unknown  colored  noise.  The  performance  of  those  algorithms  is 
frequently  very  poor  when  used  in  practical  applications.  To  derive  more  robust  al- 
gorithms, we  may  need  to  relax  the  spatially  white  noise  assumption  often  used  by 
these  algorithms.  A  robust  RELAX  algorithm  has  been  proposed  in  Chapter  3  for 
the  multi-snapshot  case  of  angle  and  waveform  estimation  of  narrowband  plane  waves 
arriving  at  a  uniform  linear  array.  The  RELAX  algorithm  does  not  assume  that  the 
noise  is  spatially  white.  The  RELAX  algorithm  is  conceptually  and  computationally 
simple;  its  implementation  mainly  requires  a  sequence  of  fast  Fourier  transforms. 
The  performance  of  the  RELAX  algorithm  has  been  evaluated  and  compared  with 
those  of  other  well-known  algorithms  including  ESPRIT  with  forward/backward  spa- 
tial smoothing,  MODE/WSF,  and  AP/ANPA  by  the  results  of  both  simulated  and 
experimental  data.  We  have  shown  that  better  signal  parameter  estimates  can  be 
obtained  by  using  RELAX  as  compared  to  using  these  other  algorithms. 

Since  most  high  resolution  narrowband  algorithms  are  not  readily  applicable  to 
wideband  sources,  focusing  matrices  are  often  needed  for  wideband  sources  and  their 
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use  can  result  in  biases  in  the  parameter  estimates.  In  chapter  4  we  have  presented 
the  wideband  RELAX  algorithm  for  the  case  of  wideband  sources.  The  wideband 
RELAX  naturally  focuses  the  narrowband  components  in  the  frequency  domain.  We 
have  used  both  numerical  and  experimental  examples  to  demonstrate  the  performance 
of  the  wideband  RELAX  algorithm  and  compared  its  performance  with  that  of  the 
CSM-ESPRIT.  We  have  shown  that  the  wideband  RELAX  algorithm  can  outperform 
the  CSM-ESPRIT  algorithm. 

To  exploit  the  advantages  of  array  geometries  and  receiving  properties  of  antenna 
sensors,  we  have  proposed  a  Co-centered  Orthogonal  Loop  and  Dipole  (COLD)  array 
in  Chapter  5.  The  Cramer-Rao  bound  for  the  COLD  array  has  been  compared 
with  that  for  the  crossed  dipole  array.  We  have  shown  that  with  the  COLD  array, 
both  angle  and  polarization  estimation  performance  can  be  significantly  improved 
as  compared  to  using  a  similar  crossed  dipole  array.  We  have  also  described  how 
to  use  the  asymptotically  statistically  efficient  MODE  algorithm  for  both  angle  and 
polarization  estimation  with  the  COLD  array.  Numerical  examples  have  been  given  to 
show  that  MODE  gives  more  accurate  angle  and  polarization  estimates  than  MUSIC 
and  NSF,  especially  for  highly  correlated  or  coherent  signals. 

A  RELAX  algorithm  has  been  devised  in  Chapter  6  for  angle,  polarization, 
and  waveform  estimation  of  narrowband  plane  waves  arriving  at  a  uniform  linear 
COLD  array.  This  RELAX  algorithm  is  also  both  conceptually  and  computationally 
simple  and  is  robust  against  the  violation  of  the  spatially  white  noise  assumption.  We 
have  presented  numerical  examples  to  demonstrate  the  performance  of  the  RELAX 
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algorithm  for  angle,  polarization,  and  waveform  estimation  with  the  COLD  array  and 
compared  it  with  the  performance  of  the  MODE  algorithm.  We  have  shown  that  this 
RELAX  can  outperform  MODE  and  yet  be  computationally  faster  than  MODE. 

7.2     Contributions 

The  main  contributions  made  in  this  dissertation  are: 

•  Developed  the  RELAX  algorithm  for  angle  and  waveform  estimation  of  narrow- 
band plane  waves  arriving  at  a  uniform  Hnear  array.  The  RELAX  algorithm  is 
robust  against  the  presence  of  unknown  spatially  colored  noise,  which  occurs 
frequently  in  practical  situations. 

•  Evaluated  the  performance  of  the  narrowband  RELAX  algorithm  and  com- 
pared it  with  those  of  other  well-known  algorithms  by  using  both  simulated 
and  experimental  data. 

•  Proposed  the  wideband  RELAX  algorithm  for  angle  estimation  of  wideband 
plane  waves  arriving  at  a  uniform  Hnear  array.  The  wideband  RELAX  does 
not  need  initial  estimates  as  other  wideband  algorithms  require  and  naturally 
focuses  the  narrowband  components  in  the  spatial  frequency  domain.  Thus  the 
wideband  RELAX  can  avoid  the  bias  problem  encountered  by  the  well-known 
coherent  subspace  methods  that  use  focusing  matrices. 
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•  Proposed  a  co-centered  orthogonal  loop  and  dipole  (COLD)  antenna  array  for 
the  estimation  of  the  parameters  of  completely  polarized  electromagnetic  nar- 
rowband plane  waves.  Due  to  the  array  structure  and  superior  receiving  pat- 
tern of  COLD,  the  performance  of  using  the  COLD  array  is  greatly  improved 
as  compared  to  using  a  conventional  crossed  dipole  array. 

•  Devised  an  efficient  MODE  algorithm  for  both  angle  and  polarization  estima- 
tion of  correlated  (including  coherent)  or  uncorrelated  incident  signals  with  a 
COLD  array. 

•  Devised  a  RELAX  algorithm  that  can  be  used  with  a  uniform  linear  COLD 
array  for  angle,  polarization  and  waveform  estimation  of  narrowband  plane 
waves. 

7.3     Future  Work 

The  topic  of  polarization  property  of  antenna  sensors  and  electromagnetic  waves 
and  polarization  estimation  has  been  discussed  by  several  researchers  [10,  31,  ?,  85, 
88].  Our  work  presented  in  Chapters  5  and  6  [34,  37]  is  the  most  recent  work  on  the 
topic  of  polarization  estimation.  All  of  these  researches  considered  only  completely 
polarized  EM  waves.  A  more  general  type  of  EM  wave  is  a  partially  polarized  EM 
wave,  whose  polarization  state  is  a  function  of  time.  Partially  polarized  waves  may 
be  found  in  many  applications  such  as  radar  and  ionospheric  radio  [41,  42].  A  large 
sample  ML  parameter  estimator  of  partially  polarized  waves  with  a  uniform  linear 
array  crossed  dipole  was  presented  in  [93].  However,  this  estimator  has  the  inherent 
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shortcoming  of  the  crossed  dipole  array,  i.e.,  the  performance  of  the  estimator  is 
sensitive  to  AOAs.  It  is  a  great  advantage  to  consider  the  parameter  estimation 
of  partially  polarized  waves  with  a  COLD  array.  Certainly,  it  is  a  interesting  and 
promising  topic  for  future  research. 

The  existing  estimation  algorithms  for  wideband  signals  have  never  been  con- 
sidered with  a  polarization  sensitive  array.  The  polarization  estimation  of  wideband 
electromagnetic  plane  waves  is  a  very  challenging  topic.  However,  it  is  possible  to 
consider  the  wideband  RELAX  algorithm  with  a  uniform  COLD  array  for  the  angle 
and  polarization  of  wideband  plane  waves.  We  shall  leave  it  as  yet  another  future 
topic. 


APPENDIX  A 
THE  CRAMER-RAO  BOUND  FOR  PARAMETER  ESTIMATES 


We  provide  herein  a  sketched  derivation  of  the  CRB  for  the  angle  and  wave- 
form estimates  of  the  incident  signals  when  the  additive  noise  is  a  zero-mean  colored 
Gaussian  random  process  with  an  unknown  covariance  matrix  Q.  The  variables  of 
the  likelihood  function  of  y(n)  include  the  unknown  elements  of  Q,  the  incident  an- 
gles, and  the  real  and  imaginary  parts  of  the  signal  waveforms.  Using  the  results  in 
[45],  we  note  that  the  CRB  of  the  signal  parameters  for  the  colored  noise  case  can  be 
obtained  from  the  CRB  for  the  white  noise  case  by  simply  replacing  A  with  Q~^/^A. 
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Also,  let 


r  =  2X;Re{P„^D^Q-^DP„}, 


(A.4) 


n=l 


where  Re{X}  denotes  the  real  part  of  X.  Define 
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where  Im{X}  denotes  the  imaginary  part  of  X.  Then  the  CR-bound  for  rj  is 
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We  remark  that  when  Q  =  a^I,  the  results  above  are  identical  to  those  in  [28]. 


APPENDIX  B 
PROOF  OF  EQUATION  (5.27) 


Let  B  be  an  I  X  (L  -  A')  matrix  whose  columns  span  the  null  space  of  A^. 
Also  let  a  (21  -  K)  x  2L  matrix  A^  be  defined  as 
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B^®I 
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We  first  show  that  the  columns  of  the  matrix  A  span  the  null  space  of  A    . 
Since  B^A  =  0  and  vf  Ujt  =  0  ,  we  obtain 


A^A 


(b^a)  ®  i)  U 


=  0. 
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This  result  shows  that  the  columns  of  A  belong  to  Af{A    ).  Moreover, 
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Thus  A  has  full  column  rank  2{L  -  K)  +  K  ^  2L  -  K.  Since 


dim 


'AfiA"")]  =  2L  -  K, 


(B.4) 


-tH. 


it  follows  that  the  columns  of  A  span  the  entire  null  space  J\f{K    ).  Hence 


H  x\-^  \H 


P^    =    Pa  =  A(A^A)"  A 

=    (B(8)I)[(B^B)"'®I](B^®I) 

+(At^  ®  I)V{V^[(A^A)"'  ®  I]v}"V^(At  ®  I).  (B.5) 


Since  Pb®i  =  Pa®I'  ^^  conclude  that  Equation  (5.27)  must  hold  true. 
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